PMF不是画出来就算数：从收敛、重加权到2D自由能面的物理判据

PMF不是画出来就算数：从收敛、重加权到2D自由能面的物理判据很多人第一次做 PMF 时，最容易掉进一个坑：图是画出来了，但物理上并不一定成立。问题在于，能画出来，和能不能当成平衡自由能解释，是两回事。这篇文章只回答几个更基础、也更容易出错的问题：已有数据什么时候足够支持 PMF，什么时候只能报局部结果，什么时候必须重加权，什么时候 2D 图虽然能画，但其实不该把它写成“收敛的自由能面”。结论 PMF 的定义本身并不难，真正困难的是采样是否真的支持这个定义。无偏 MD 确实可以直接给自由能，但前提是分析段已经平稳，而且目标坐标空间被充分访问；只要存在偏置、约束、umbrella 或多窗口合并，就不能跳过重加权。 2D PMF 不是“多画一个维度”那么简单，而是对采样混合提出了更高要求。如果某些区域从来没被访问过，任何后处理都不能把真实自由能补出来；因此，很多时候你真正能安全报告的，并不是全局 PMF，而是局部 PMF、条件分布或状态占据。 PMF 到底是什么对一个集合变量 $\xi$，平衡自由能剖面定义为： \[F(\xi) = -k_B T \ln P(\xi) + C\] 如果有两个集合变量 $\xi,\eta$，对应的二维自由能面就是： \[F(\xi,\eta) = -k_B T \ln P(\xi,\eta) + C\] 公式的通俗解释这两个式子真正表达的是一句很朴素的话：某个状态如果在平衡系综里更常出现，它的自由能就更低。所以，问题的核心从来不是“会不会取负对数”，而是你算出来的 $P(\xi)$ 或 $P(\xi,\eta)$ 到底是不是平衡分布，这个分布覆盖的是全局空间还是只覆盖了一个局部盆地，以及每个 bin 里到底有多少有效独立样本。这三件事，才真正决定了你的 PMF 能不能被当成物理结果来解释。在后面的例子里，我会经常用 P2 和 Z 这两个符号。这里可以先把它们通俗地理解成两类常见坐标：P2 代表某种取向序参量，也就是“分子更偏向平躺、倾斜还是竖直”的量化描述；Z 代表某种位置坐标，例如分子相对于界面、膜中心或参考平面的距离。你完全可以把它们替换成自己体系里真正关心的两个集合变量。什么叫“物理上正确”的 PMF 如果想让一条 PMF 在文章里站得住脚，至少要同时满足四件事：数据来自同一个目标系综用来分析的轨迹段已经进入平稳区你关心的坐标范围内发生了足够的往返跃迁误差估计使用的是有效样本数，不是总帧数只要这四条里缺一条，图可能仍然能画出来，但解释时就必须明显降级。第一关：是不是同一个统计系综这一点最容易被忽视。如果所有数据都来自同一统计系综，也就是温度一致、压力设置一致、力场和拓扑一致、体系组成与边界条件一致，同时没有额外偏置或约束，那么这些轨迹才有资格被当作同一个平衡分布的样本来合并分析。那么你可以直接从直方图或核密度估计（KDE）得到 $P(\xi)$，再转成自由能。但只要出现下面任一种情况，就不能把所有帧直接混在一起做直方图：情况为什么不能直接混合对某个坐标加了 umbrella 势采样分布已经被显式改权，不再对应原始无偏分布加了位置约束或取向约束体系访问相空间的方式被限制，直方图不再代表自然占据做过 steered MD 或 pulling 轨迹带有外场驱动，不能直接当成平衡样本合并了不同温度的数据不同温度对应不同平衡分布，不能简单拼接合并了不同哈密顿量或不同参数的数据势能面本身不同，统计权重自然也不同这时你要处理的已经不是“无偏概率”，而是“被改权重后的采样概率”。必须重加权，常见工具就是 WHAM、MBAR，或者更一般的重加权流程。第二关：轨迹是不是已经进入平稳区很多 PMF 最大的问题，不是采样短，而是前半段根本还没平衡。比如系统一开始从某个强行构建的初始构型出发，前几十纳秒甚至更久都还在弛豫。如果把这一段直接并进统计，得到的就不是平衡分布，而是“初始条件残留 + 平衡波动”的混合物。一个实用做法，是先做平衡段检测，再决定从哪里开始统计。常用工具是 pymbar.timeseries。这里输入的数据，不是什么特殊格式文件，而是某个集合变量随时间变化的一列数据，最常见的就是 P2(t) 或 Z(t) 这样的时间序列： python - <<'PY' from pymbar import timeseries import numpy as np P2_t = np.loadtxt('P2_t.dat') t0, g, Neff = timeseries.detect_equilibration(P2_t, nskip=10) print(t0, g, Neff) PY 如果你手里保存的是多列文件，例如同一份文件里同时有时间、P2 和 Z，那就应该先把你想分析的那一列取出来，再送进 detect_equilibration()，而不是把整张表不加区分地直接读进去。这里最值得报告的，不是“我跑了多少 ns”，而是平衡起点 $t_0$、统计低效因子 $g$ 和有效样本数 $N_{\mathrm{eff}}$。真正决定误差条大小的，是独立样本有多少，不是帧有多少。很多时候看起来“已经有几十万帧”，但如果自相关很强，真正能用于统计判断的独立样本可能并不多。第三关：有没有真正发生“来回走动” 这是判断 PMF 是否可信的核心。真正有用的判断，不是“分布看起来挺宽”，而是体系有没有在你关心的几个主要状态之间真正来回走动，也就是是否发生了足够多的往返跃迁（round trips）。对 1D 和 2D PMF，要求到底差在哪里目标至少要看到什么不能轻易下的结论 1D PMF 主要盆地被多次访问，盆地之间有往返跃迁，不同重复给出相近边缘分布只有单盆地波动时，不应宣称得到全局 PMF 2D PMF 两个坐标都被实质性访问，且在固定第一维时第二维也能混合，不同区域之间整体连通如果第二维几乎没动，或固定某一维后另一维几乎不跨峰，就不应宣称得到全局 2D 自由能面如果体系只在一个盆地附近晃动，那么你当然也能画出一条曲线，但那更接近“局部热涨落的自由能近似”，而不是全局 PMF。二维情况则更严格，因为它要求你不仅采到 $\xi$，还要在不同 $\eta$ 条件下把 $\xi$ 也采匀；一旦第二维只是窄范围波动，这张 2D 图通常就只能算局部地形。一个最常见的误区：能画 2D，不等于应该发 2D 很多人会这样做：选两个坐标，做二维直方图，再对联合概率取负对数，最后得到一张彩色图。从程序角度看完全没问题，但从物理角度看，可能只说明一件事：你的轨迹在一个局部区域里留下了很多点。这时真正应该问的，不是“图是不是好看”，而是三个更扎实的问题。第一，第二维是不是只覆盖了一个很窄的范围；如果是，那么 2D 图只是把局部波动展开成二维，并没有真正回答更大的自由能问题。第二，高自由能区域是“真的高”，还是“根本没采到”；没有访问到的格点，在视觉上很容易被误读成高能区，但统计学上它可能只是空白区。第三，盆地之间的通道是物理能垒，还是统计断裂；如果两个盆地中间几乎没有过渡点，你看到的未必是高能屏障，也可能只是采样没有连通，更专业地说，就是这些区域之间缺少足够的统计连通性。如果这些问题答不上来，最稳妥的表述通常不是“得到了全局 2D PMF”，而是把口径主动降到“局部 2D 自由能地形”“条件分布 $P(\xi\mid\eta)$”或者“已结合区间内的取向自由能”。什么时候无偏 MD 足够无偏 MD 适合回答的问题，其实比很多人想象得更有限，但也更扎实。与其笼统地说“能不能算 PMF”，不如先区分你到底想回答哪一类问题。目标无偏 MD 的适用性更合适的表述单个坐标的 1D 边缘自由能较好 1D PMF 某个局部区域内的自由能起伏较好局部 PMF 分箱后的状态占据比较较好条件分布或占据统计跨多个盆地的全局自由能谨慎只有在多次跨盆地跃迁后才可报告同时含位置与取向的 2D 自由能面很谨慎通常先降级为局部 2D 或条件分布含解离、再结合、重排等慢过程很谨慎往往需要增强采样支撑如果你的无偏轨迹从头到尾都没有离开某个状态盆地，那么最合理的结论不是“体系没有别的态”，而是：当前采样没有能力回答这个问题。什么时候必须用 WHAM 或 MBAR 这个判断其实很干脆：只要采样权重被改过，就要重加权。与其把这一条说成一句口号，不如直接看常见场景：场景能不能直接做直方图推荐处理同一无偏 MD 可以直方图或 KDE umbrella 窗口不可以 WHAM 或 MBAR 多温度数据合并不可以 MBAR 有约束或 pulling 不可以显式重加权多个偏置窗口做 2D 分布不可以先去偏，再做联合分布如果你手上已有沿某个坐标布置好的 umbrella 窗口，那么它们通常足够支持可靠的 1D PMF。至于能不能进一步得到 2D PMF，要看另一个坐标在每个窗口里是不是也混合得足够好。主坐标被偏置采到，并不自动意味着旁观变量也已经收敛，这一点在实际分析里经常被误判。一个非常实用的判断：你到底能安全声称什么诊断结果最稳妥的说法只有一个局部盆地被采到局部自由能或局部涨落 1D 有多次跨峰跃迁，重复一致可以报告 1D PMF 2D 中第二维很窄只报告条件分布或局部 2D 地形 umbrella 在主坐标重叠良好，但副坐标混合差主坐标 PMF 可信，2D 结果仅作定性参考每个窗口内副坐标多次跨峰，重复一致可以认真讨论 2D PMF 这张表背后的原则其实很简单：结论的口径，必须和采样能力匹配。很多结果并不是“完全不能发”，而是应该主动把口径降到“局部 PMF”“条件分布”或者“占据统计”这一层，这样反而更稳。收敛不能只看“曲线变平” 很多人判断收敛时，只看 PMF 曲线后半段是不是“不怎么变了”。这远远不够，因为一条表面平滑的曲线，可能只是建立在高度相关、重复不一致、或者根本没有跨盆地跃迁的数据上。更可靠的收敛证据链更可靠的判断，通常要把下面几类证据合在一起看：先看结果会不会随时间继续漂，也就是是否仍在发生系统性漂移；再看不同重复是否支持同一组物理结论；接着看你到底有多少真正独立的样本；最后再确认主要状态之间有没有真正发生来回切换，也就是是否存在足够的往返跃迁。时间分块分析：把前 1/3、前 2/3 和全部数据分别算一次 PMF。这样做的目的，不是为了多画几条线，而是看结果会不会继续变。如果主要盆地位置、相对深度和势垒高度还在系统性漂移，那就说明体系还在持续演化、尚未真正稳定下来，此时“看起来平滑”并不等于已经收敛。重复一致性：不同重复轨迹给出的分布或 PMF 应该大体一致。这里最重要的不是三条线能不能完全重合，而是它们是否支持同一个物理结论。如果不同重复之间差异明显，最常见的解释不是“体系本来就这样”，而是混合仍然不足，也就是每条轨迹还在各自记着不同的初始路径。自相关分析：报告 $g$ 和 $N_{\mathrm{eff}}$，确认自己不是在用几十万帧去假装拥有几十万个独立样本。连续轨迹里的相邻帧往往很像，所以“帧数很多”不等于“信息很多”。这一步本质上是在修正相关样本导致的误差低估，也就是给误差条去水分，说明到底有多少真正能独立贡献统计信息的数据点。跃迁计数：主要盆地之间要有实质性的往返，而不是只在一个盆地里高频抖动。很多人看到时间序列很活跃，就以为体系采样得很好，但如果这些波动始终发生在同一个局部盆地里，那么关键状态之间的相对自由能差其实还没有被真正比较过。没有跨盆地跃迁时，很多相对自由能差并不稳。窗口重叠：对 umbrella 来说，相邻窗口必须足够连通。如果相邻窗口之间几乎没有共同覆盖的区域，WHAM 或 MBAR 就很难把整条 PMF 稳稳地拼起来。这时数学上虽然还能算，物理上却可能只是把几段彼此脱节的局部结果硬接在一起；更规范地说，就是窗口之间缺少足够的概率分布重叠。 umbrella 数据至少要看什么对于 umbrella，gmx wham 的常规检查项很重要： gmx wham -it tpr-files.dat -if pullf-files.dat -o pmf.xvg -hist hist.xvg -ac 这里至少要看三件事，而且最好把它们理解成“这条 PMF 能不能被顺畅接起来”的三个层次检查：相邻窗口直方图有没有足够重叠。这是最基础的一关。如果相邻窗口几乎不相交，那么后处理再漂亮，也只是把统计上彼此脱节的区间强行缝在一起，整条曲线会缺少真正的连接。自相关时间是不是已经大到接近单窗口长度。这一步是在问：单个窗口里到底有没有采到足够多的独立信息。如果一个窗口里有效独立样本本来就很少，那么它对整条 PMF 的贡献会既不稳定又很难估误差；此时窗口数量再多，也不等于每个窗口都真的达到局部统计稳定。不同窗口拼起来后有没有明显断链。所谓断链，不一定表现成肉眼可见的大跳跃，也可能表现为某些区间误差异常、重复不一致，或者对分析参数极其敏感。如果一条 PMF 只要稍微改一下 bin、平滑或截断方式就明显变样，那通常不是“图画风不同”，而是底层采样还不够扎实。如果某些窗口几乎没有重叠，或者窗口内采样时间和自相关时间是一个量级，那这套 PMF 就很难让人放心。 2D PMF 什么时候才值得做更关键的问题是：什么时候做 2D PMF 比做 1D 或条件分布更有信息增益。通常至少要同时满足三点：两个坐标都对应你真正关心的慢过程，这两个坐标在数据里都被实质性采样到了，而且在固定第一维时第二维不是“卡死”的，也就是没有被困在某个狭窄取值范围里。少了其中任何一条，二维分析带来的往往不是新信息，而是新噪声。如果不满足，2D 往往只会带来两个后果：图更花哨，误差更大。因为二维一上来就会遭遇“维数灾难”：格点数一多，平均到每个 bin 的有效样本数会迅速下降，空 bin 和噪声会明显增加。所以，在下面这些情况下，不做 2D 反而更专业：如果第二维只是辅助解释变量，如果第二维的采样范围很窄，如果第二维的混合时间明显比单窗口长度更长，或者你的核心结论本质上靠 1D 就已经成立，那么继续硬做 2D 往往只会增加图的复杂度，而不会提高结论的可信度。还有一个细节：有些序参量自带“几何熵” 如果你用的是角度、取向序参量，或者由角度变换得到的量，那么要小心一个问题：原始分布里可能混进了变量测度本身带来的偏置。最直观的例子就是方向相关变量。即使体系完全各向同性，某些取向序参量的概率分布也未必是均匀的。这意味着直接计算 \[F(\xi) = -k_B T \ln P(\xi) + C\] 得到的可能既包含真实相互作用偏好，也包含“随机几何本来就更容易落在某些值附近”的贡献。这时最常见的处理方式有两种：报告方式含义适合的讨论场景原始 PMF 包含变量测度带来的几何熵讨论状态占据、总体分布相对参考分布的超额自由能更突出相互作用导致的偏好讨论取向偏好、界面诱导效应这不是所有体系都必须做，但如果你的核心结论高度依赖“取向偏好”，那这个问题最好提前想清楚。否则读者看到的“最低谷”，有一部分可能只是变量定义自带的几何效应，而不全是体系相互作用本身。一个面向实战的工作流 graph TB A["拿到已有轨迹"] --> B["先分清：无偏数据还是有偏数据"] B --> C["确定目标：1D、局部2D、还是全局2D"] C --> D["检测平衡段：t0、g、Neff"] D --> E["检查跃迁、重复一致性、窗口重叠"] E --> F{"采样是否支持目标结论"} F -->|支持| G["报告 PMF，并给出误差与收敛证据"] F -->|部分支持| H["降级为局部 PMF、条件分布或状态占据"] F -->|不支持| I["补采样或重新设计增强采样方案"] 这个流程最重要的一步，不是“画图”，而是中间那个判断：采样能力到底支不支持你想说的话。真正成熟的分析，不是把所有图都画出来，而是知道哪些图值得认真解释，哪些图只能当辅助材料。结果该怎么讲，才更站得住脚一张自由能图要站得住脚，关键不在于修饰，而在于先把哪里可信、哪里还不能多说讲清楚：先说明平衡段和有效样本是怎么处理的。如果一开始就交代你已经剔除了前期非平衡部分，并且按相关性修正了有效样本数，读者会更容易接受后面的自由能结果，因为他知道这些曲线不是把所有帧不加区分地堆出来的。再说明 1D 结果为什么可信。如果主要状态之间已经出现多次往返跃迁，而且不同重复支持同一个结论，那么这时去讨论 1D PMF 的相对高低才更有底气，因为它背后有明确的动力学采样证据。谈到 2D 结果时主动限定范围。如果二维图只有一部分区域采样得比较扎实，那就只讨论那一部分，把它明确写成局部自由能地形或条件分布。这样做不会削弱文章，反而会让读者觉得你的判断更稳。对空白区和混合不足区保持克制。没有访问到的区域就不要硬解释，混合明显不足的方向也不要勉强下定量结论。这样做不是示弱，而是在保护结论的可信度。这种写法的价值不在于“更谨慎”，而在于把真正确定的部分讲扎实，把暂时不能确定的部分老老实实留白。最后总结 PMF 真正难的地方，从来不是软件命令，而是你是否对“这张图能回答什么问题”有清醒判断。无偏 MD 确实可以直接给自由能，但前提是轨迹分析段已经平稳、混合、可重复。如果连主要状态之间的往返都没有发生，那么图上看到的更多只是局部波动，而不是可以放心解释的全局自由能。只要数据里存在偏置、约束、umbrella 或多窗口拼接，就必须认真做重加权。这不是后处理里的可选美化步骤，而是把“被改过权重的采样”还原成目标分布所必需的物理操作。 2D PMF 的门槛显著高于 1D PMF，因为它要求两个坐标都被充分访问，而且在固定其中一维时另一维也要发生足够混合。很多 1D 看起来已经稳定的数据，一到二维分析就会暴露出空白区、断裂区和高噪声问题。没采到就是没采到，后处理不能替代真实采样。无论是更平滑的直方图、更复杂的重加权，还是更漂亮的二维彩图，都不能凭空恢复从未被访问过的状态或通道。当采样只支持局部结论时，老老实实报告局部结论，反而更有说服力。把结果写成局部 PMF、条件分布或状态占据，通常比强行宣称“全局自由能面已经收敛”更专业，也更经得起追问。如果把这套判断标准先建立起来，你之后无论做无偏 MD、umbrella、metadynamics，还是更复杂的多维自由能分析，很多技术决策都会清楚得多。

Molecular Dynamics · 2026-03-31

BioEmu能把蛋白动力学采样推多远：激酶成功，转运体与隐蔽口袋暴露边界

BioEmu能把蛋白动力学采样推多远：激酶成功，转运体与隐蔽口袋暴露边界本文信息标题：Accelerated sampling of protein dynamics using BioEmu augmented molecular simulation 作者：Soumendranath Bhakat，Eva-Maria Strauch 发表时间：2026年2月21日（bioRxiv 预印本）单位：AlloTec Bio Inc.（美国密苏里州圣路易斯）；Washington University in St. Louis School of Medicine, Division of Infectious Diseases（美国密苏里州圣路易斯）引用格式：Bhakat, S., & Strauch, E.-M. (2026). Accelerated sampling of protein dynamics using BioEmu augmented molecular simulation. bioRxiv. https://doi.org/10.64898/2026.01.07.698041 源代码与相关工具： BioEmu：https://github.com/microsoft/bioemu H-packer：https://github.com/gvisani/hpacker CryoPhold：https://github.com/strauchlab/cryoPhold MDML：https://github.com/svats73/mdml/tree/main 摘要这篇预印本提出了一条把生成式AI构象生成、无偏分子动力学模拟和马尔可夫状态模型串起来的工作流。作者先用 BioEmu 生成蛋白质骨架构象，再补全侧链、做慢特征分析与聚类，最后从代表性结构出发跑多条短程 MD，并用 MSM 恢复符合玻尔兹曼权重的构象分布。在 CDK2 与 BRAF 这类丝氨酸／苏氨酸激酶上，这条路线确实能捕获 DFG-in 到 DFG-out 的稀有转变，还能解析 V600E 突变诱导的群体迁移。更进一步，作者把 BioEmu 与 Cryo-EM 重加权结合，用于构建 GlyT1 的全原子构象系综。不过，论文同样强调了一点：BioEmu 并不是普适的动力学万能钥匙。在 GlyT1 与 PlmII 这类强依赖侧链构象异质性的体系里，BioEmu 派生的初始系综并没有覆盖足够广的功能相关状态，后续 MD 也就难以“凭空补回来”。核心结论 BioEmu 加短程 MD在激酶体系里确实有效，能用累计 5 μs 的模拟捕获 DFG-in 到 DFG-out 转变，而对照的 rMSA-AF2 路线即使做到 8 μs 仍主要困在 DFG-in 这套方法不只是找到“终态”，还能够解析中间态、亚态和群体比例，例如 CDK2 激活环折叠／伸展状态与 BRAF 的 DFG-Phe 旋转异构体分布对 V600E BRAF，方法成功恢复了突变诱导的群体转移，包括 DFG-Phe 从 PheF1 向 PheN 的偏移，以及 αC 螺旋向更活性样构象偏移把 BioEmu 与 Cryo-EM 贝叶斯重加权结合后，可以得到 GlyT1 的全原子先验系综，但采样仍然不完整，尤其是 inward 态与 Y62 翻转论文最重要的结论其实是边界条件：当动力学高度依赖侧链异质性时，只有骨架多样性往往不够，BioEmu v1.0 的优势会明显下降背景蛋白质功能往往不是由单一静态结构决定的，而是由多个亚稳态之间的相对群体与相互转化共同决定。对药物研发来说，这一点尤其关键，因为变构口袋开放、激活环重排、跨膜转运开关、蛋白—蛋白相互作用界面暴露，很多都属于低概率但功能关键的稀有事件。传统无偏 MD 最大的问题是时间尺度。很多功能相关转变隔着很高的自由能垒，常规模拟在可接受的算力预算内根本跨不过去。增强采样方法当然能帮忙，但常常要提前指定集体变量，或者引入偏置势，后续还得重新加权。问题不在于这些方法不好，而在于它们通常依赖较强的经验判断。这几年生成式 AI 进入分子模拟领域后，一个自然的问题是：能不能让 AI 先把构象空间“撒开”，再由物理模拟去恢复真实分布？这篇文章的思路正是如此。不过作者没有把 BioEmu 包装成万能替代品，而是很认真地比较了它在不同体系中的表现，最后给出的结论是：它在某些问题上很强，但也有非常具体、非常物理的失效场景。关键科学问题 BioEmu 生成的构象系综，能不能真正作为稀有动力学事件的高质量初始分布把 BioEmu、短程 MD 与 MSM 串起来之后，能否恢复有物理意义的平衡态群体与自由能面这套路线在成功体系与失败体系之间，分界线到底在哪里如果体系的关键转变高度依赖侧链翻转、局部闸门残基或隐蔽口袋开启，BioEmu 是否还足够好用创新点提出了一个相对清晰的两阶段流程：先用 BioEmu 做广覆盖，再用 MD＋MSM 做物理校正不只展示成功案例，还专门纳入 GlyT1 和 PlmII 这类有挑战的反例体系把 CryoPhold 的 Cryo-EM 贝叶斯重加权流程扩展到 BioEmu 先验系综用同一篇文章同时回答“这方法什么时候有效”和“什么时候会失灵”，这一点其实比单纯展示漂亮案例更有价值研究内容整体工作流：先铺开构象，再交给物理学筛选图1：BioEmu 种子分子模拟的整体工作流输入是蛋白质序列，BioEmu 先生成约 500 个仅含骨架的单体构象 H-packer 负责补全侧链，把骨架系综转换成全原子表示作者对 Cα–Cα 距离做慢特征分析，并在前两个慢特征上进行 K-means 聚类，得到 50 个代表性结构这 50 个代表性结构分别启动 100 ns 无偏 MD，总计 5 μs 所有轨迹最后交给 MSM 统一整合，输出自由能面、宏观态群体和亚态分布这张图的重点不是“AI 替代了 MD”，而是AI 改变了初始结构分布。作者反复强调，后面的 MD 仍然是物理驱动的，只不过 BioEmu 提供了一个更可能覆盖稀有态的起点。方法的关键逻辑：覆盖率先行，但物理意义不能省这篇文章最值得记住的一句话可以概括成：BioEmu 负责把你带到更多地方，MSM 负责告诉你哪些地方真正重要。如果只看 BioEmu 本身，它给出的是构象多样性，而不是严格的平衡分布。作者因此没有直接把 BioEmu 输出当答案，而是把它当作更聪明的初始构象提案器。后续的全原子 MD 和 MSM，才是赋予这些结构统计物理意义的步骤。这一点也解释了为什么作者坚持用对照组。文章不是简单展示“BioEmu 能采到什么”，而是要比较：同样是短程无偏 MD，不同初始构象覆盖到底能把结果拉开多大差距。激酶测试：BioEmu 的最佳表现出现在 DFG 翻转问题上图2：MSM 加权自由能面解析 BRAF 与 CDK2 的 DFG-in 到 DFG-out 转变 A、C 是 BioEmu 种子模拟得到的自由能面，分别对应 apo BRAF 与 apo CDK2 B、D 是 rMSA-AF2 增强 MD 的对照结果黑点是初始构象系综投影，作者用它来直观看出初始覆盖范围 E 给出了 DFG-in 与 DFG-out 的代表性结构，salmon 色对应 DFG-in，cyan 色对应 DFG-out，重点看的是 DFG-Phe、Lys、Glu 的相对位置变化这组结果非常直观。BioEmu 种子模拟不只是跑出了更散的点云，而是真正在自由能面上覆盖到了从 DFG-in 到 DFG-out 的过渡区域。相比之下，rMSA-AF2 的初始系综和后续模拟几乎都局限在 DFG-in 附近。这里最有说服力的不是“总能量更低”之类抽象说法，而是一个非常实际的比较：BioEmu 路线总模拟时间是 5 μs，对照路线是 8 μs，但后者仍没能真正跨出 DFG-in 盆地。这说明在这类问题上，初始构象覆盖确实比单纯延长短程模拟更重要。 CDK2：不仅采到 DFG-out，还采到了更细的活化相关异质性图3：BioEmu 增强模拟解析 apo CDK2 的 DFG-Phe、αC 螺旋与激活环亚态 A 是 DFG-in 宏观态内不同 DFG-Phe 旋转异构体，以及 αC 螺旋 LGL／LGU 和激活环 ACin／ACout 的相对群体 B 把激活环距离投影到 DFG 相关的两个距离坐标上，显示 DFG-out 更偏向折叠激活环 C 叠合了代表性 DFG-in 与 DFG-out 结构，突出显示DFG-Phe 翻转与激活环折叠如果图2告诉读者“BioEmu 能跨盆地”，那图3告诉读者的是：它不只会跨盆地，还能把盆地里的精细异质性解析出来。在 apo CDK2 里，作者不仅看到了 DFG-in 与 DFG-out 两个终态，还看到了 DFG-in 内部的不同 DFG-Phe 亚态，以及 αC 螺旋与激活环的耦合变化。尤其是从 DFG-in 到 DFG-out 时，激活环从 ACout 向 ACin 转移，这正是从更活性样构象走向更非活性样构象的重要标志。换句话说，BioEmu 的价值不只是“帮忙见到稀有终态”，而是能让后续 MSM 在更合理的初始覆盖上，恢复出与功能转换相关的层级化构象景观。 V600E BRAF：群体转移而不是单一结构切换，才是更难也更有用的测试图4：V600E 突变如何把 BRAF 系综推向更活性样构象左侧柱状图比较野生型与 V600E 在 DFG-in 宏观态内的 PheN、PheF1、PheF2 群体中间柱状图比较 αC 螺旋在 LGL 与 LGU 两种构象下的群体变化右侧结构示意图标出 Phe595、Lys483、Glu501，并用蓝色与米色展示更偏 DFG-in／DFG-out 或 LGL／LGU 的构象差异这一部分是全文最接近“生物学解释”的地方。作者并不是简单说 V600E 更活跃，而是用群体分布具体展示：在 DFG-in 宏观态内部，V600E 会把 DFG-Phe 的侧链旋转异构体从 PheF1 推向 PheN。同时，αC 螺旋也更偏向“in”状态，也就是 LGL。这很重要，因为突变激活常常不是把蛋白从一个完全静止的构象“掰”到另一个，而是让整个系综在多个亚态之间重新分配权重。这篇文章的一个亮点就在于，它确实把这种“群体转移”用 MSM 权重给量化了出来，而不只是画一张构象示意图就结束。把 Cryo-EM 和 BioEmu 接起来：GlyT1 是更接近真实应用场景的测试图5：BioEmu 先验系综经 CryoPhold 贝叶斯重加权后，得到 GlyT1 的全原子构象集合左侧是原始 BioEmu 系综和 SFA 聚类后的 50 个代表性结构右上是三张 Cryo-EM 参考图，对应 inward、occluded 与 outward 三种状态，分辨率分别约为 3.35 Å、2.58 Å 和 3.22 Å 右下是重加权后的全原子 CryoPhold 系综，橙色、青绿色、紫色分别对应 inward、occluded、outward 这部分很值得关注，因为它把“AI 给先验 + 实验给约束 + MD 给动力学”这三件事真正串了起来。作者不是直接拿 BioEmu 去解释 GlyT1，而是先通过 Cryo-EM 参考图做贝叶斯重加权，得到更接近实验的全原子后验系综。从概念上说，这一步很漂亮：BioEmu 给广覆盖的起点，Cryo-EM 给状态约束，CryoPhold 把两者合成更可信的结构先验。如果只看工作流设计，这其实是全文最有方法学延展性的部分。但问题也从这里开始：GlyT1 并没有被完全采开图6：在 GlyT1 上，BioEmu 系综的覆盖不足开始暴露出来 A 标出 GlyT1 的关键热点残基，尤其是 Y62、W322、R71、D474，它们共同定义了状态转变相关的局部几何 B 是 BioEmu 种子模拟在 TM1–TM6 与 TM1–TM10 距离空间中的采样结果 C 是 rMSA-AF2 种子模拟的对照，明显覆盖到更多 inward、occluded、outward 区域 D、E 则比较了 Y62 的 χ1／χ2 二面角采样，显示 BioEmu 路线对 Y62 翻转的覆盖明显不足这张图非常关键，因为它直接告诉读者：BioEmu 并不是在所有体系里都比 rMSA-AF2 更强。在 GlyT1 中，作者发现 CryoEmu 增强模拟虽然能较好采到 outward 与 occluded，但对 inward 态以及 Y62 翻转的恢复并不充分。这个结果和前面激酶体系的成功形成鲜明对比，也说明 GlyT1 的关键动力学更依赖局部残基闸门与侧链重排，而不只是主链骨架的大尺度移动。也就是说，对某些跨膜转运体来说，单纯把骨架铺得更开并不够。真正控制状态切换的，可能是像 Y62 这样的局部“盖子”残基，而这恰恰是 BioEmu v1.0 不擅长的地方。 PlmII：隐蔽口袋开启再次证明，侧链问题绕不过去图7：在 PlmII 的隐蔽口袋开启问题上，rMSA-AF2 反而明显优于 BioEmu A 是 BioEmu 增强模拟得到的 Trp41 χ1／χ2 自由能面，基本只覆盖主态 B 是 rMSA-AF2 的对照结果，可以看到更多离散盆地，其中圈出的区域对应隐蔽口袋开启相关状态 C 给出 Trp41 翻转的结构示意，说明这个侧链运动与口袋暴露直接相关如果说 GlyT1 已经让人开始怀疑“骨架覆盖是否足够”，那 PlmII 几乎就是把这个问题钉死了。作者明确指出，PlmII 的隐蔽口袋开启依赖 Trp41 侧链翻转，而 BioEmu 生成的初始系综在这件事上的构象多样性太有限，所以后续 MD 也很难补救。这也是全文最值得记住的负面结论之一：对由关键侧链翻转主导的构象开关，BioEmu v1.0 的瓶颈不在后续采样，而在起跑线就没有把相关侧链异质性准备好。这篇文章真正回答的问题：什么时候该用 BioEmu，什么时候要谨慎综合激酶、GlyT1 和 PlmII 三类体系，这篇文章给出的不是一个简单的“好用／不好用”结论，而是一个更细的经验判断。更适合 BioEmu 的情形通常有这些特征：关键转变主要表现为骨架层面的宏观构象重排稀有态虽然难采，但可以由较广的主链分布触达后续短程 MD 加 MSM 足以把这些状态重新赋予物理权重相对不利的情形则包括：关键动力学由局部侧链翻转控制功能相关状态依赖少数残基构象的精细组合起始系综如果没有覆盖这些局部侧链模式，后续无偏 MD 很难在短时间内补齐这也是作者为什么会在摘要和讨论里都强调，BioEmu 更像是一个很强的构象覆盖工具，而不是自动恢复全部真实动力学的黑箱。方法细节：这套流程到底是怎么落地的体系生成与聚类 BioEmu v1.0 为每个体系生成约 500 个构象 H-packer 补全侧链，得到全原子结构作者对 Cα–Cα 距离做慢特征分析，并在前两个慢特征上进行 K-means 聚类，得到 50 个聚类中心 SFA 与聚类使用的是 MDML 软件包对 GlyT1，作者再把这 50 个聚类中心作为 CryoPhold 的先验，用于针对三张 Cryo-EM 图的贝叶斯重加权分子模拟参数使用 Amber2022 中的 tleap 进行体系准备蛋白力场是 AMBER ff14SB 水模型是 TIP3P 使用截角八面体水盒，蛋白到盒边界最小缓冲为 10 Å 先做受限最小化，再做全体系无约束最小化 Amber 拓扑通过 ACPYPE 转到 GROMACS 格式后续模拟在 GROMACS 2022 中进行体系从 0 K 升温到 300 K，先进行 500 ps NVT 升温，再进行 200 ps NPT 平衡生产模拟为无偏 100 ns，轨迹每 10 ps 保存一次温控采用 velocity-rescale thermostat，压强控制采用 Parrinello–Rahman barostat 非键相互作用截断为 1.0 nm，长程静电采用 PME，含氢键长通过 LINCS 约束 MSM 构建 MSM 使用 PyEMMA 构建激酶体系使用图2中的两个距离来区分 DFG 态 GlyT1 则使用能区分 inward、outward、occluded 的距离变量来建模从技术路线看，这篇工作的核心不在于发明了新的采样偏置算法，而在于把生成式构象先验、全原子 MD、MSM 与实验约束拼成了一条相对简洁、可复用的流程。 Q&A Q1：为什么 BioEmu 在激酶上明显成功，但在 GlyT1 和 PlmII 上表现变差？ A1：一个核心区别是，激酶 DFG 转变更多体现为主链与局部二级结构层面的构象重排，而 GlyT1 的 Y62、PlmII 的 Trp41 都属于关键侧链闸门残基。BioEmu v1.0 只显式生成骨架，侧链是后补的，所以一旦功能动力学高度依赖侧链异质性，起始覆盖就会受限。 Q2：这篇文章是在说 BioEmu 比 rMSA-AF2 更好吗？ A2：不是简单的“更好”，而是不同体系各有胜负。在 BRAF 和 CDK2 这类激酶上，BioEmu 的构象覆盖明显更广；但在 GlyT1 与 PlmII 上，rMSA-AF2 反而给出了更好的功能相关采样。作者真正想说明的是：初始系综的质量必须和问题类型匹配。 Q3：为什么作者要坚持在后面再跑 MD 和 MSM，不能直接分析 BioEmu 输出吗？ A3：因为 BioEmu 给的是结构覆盖，不是严格的平衡分布。后续 MD 才提供局部物理松弛，MSM 才负责把多条轨迹整合成有统计力学意义的群体分布与自由能面。没有这一步，BioEmu 更像“候选构象生成器”，而不是完整的动力学答案。 Q4：这项工作对药物发现最直接的启发是什么？ A4：它说明了一个很实际的策略：如果目标体系的关键动力学主要由骨架级别的大构象转变主导，BioEmu 这类模型可以显著提高稀有态触达率；但如果问题核心是局部侧链翻转、闸门残基摆动或隐蔽口袋开启，就不能指望只靠骨架多样性解决问题，必须考虑更强的侧链建模或额外实验约束。关键结论与批判性总结这篇文章最重要的价值这篇文章真正有价值的地方，不是单纯展示“AI 让采样更快”，而是把这个命题拆开讲清楚了。作者既给出了 BioEmu 在激酶体系上的漂亮成功案例，也非常诚实地展示了它在 GlyT1 和 PlmII 上的失败边界。这种写法反而更有参考意义，因为它把方法的适用前提说透了。主要优点成功案例很有说服力：BRAF 与 CDK2 的 DFG 转变确实被采到了，而且对照组差距明显不只看终态：文章分析了中间态、亚态、群体分布和突变诱导的群体转移，信息密度很高工作流具有可操作性：BioEmu、H-packer、MDML、GROMACS、PyEMMA、CryoPhold 串起来后，路线相对明确对失败模式有清楚归因：作者把问题聚焦到侧链异质性不足，这个解释既具体又有物理直觉局限性 BioEmu v1.0 不显式建模侧链，这会直接限制对 Y62、Trp41 这类关键残基翻转的覆盖当前流程主要面向单体蛋白，对蛋白—蛋白或蛋白—配体体系的适用性仍有限虽然结果与已知机制一致，但很多系统仍缺少更直接的实验定量验证成败在很大程度上取决于初始系综是否覆盖到真正相关的局部自由度，这意味着方法仍然需要系统特异性判断对后续工作的启发如果未来的生成模型能更好处理全原子级别的侧链异质性，这条路线的适用范围会明显扩大把 Cryo-EM、DEER、FRET 等实验信息与生成模型输出做更紧的耦合，可能是提高可靠性的关键方向对于隐蔽口袋和局部闸门问题，后续方法很可能需要从“只学骨架”走向同时学习骨架与关键侧链坐标总体来看，这篇文章最值得记住的 punchline 不是“BioEmu 已经解决了蛋白动力学采样”，而是：BioEmu 确实能显著改善一类问题，但它的边界恰好暴露了下一代生成式分子模拟模型最该补的课。

Molecular Dynamics · 2026-03-18

神经关系推断：从MD轨迹中学习蛋白质长程变构相互作用

Molecular Dynamics · 2026-01-25

变构激活的动态基础：恶性疟原虫蛋白激酶G的长程通信机制

Molecular Dynamics · 2026-01-22

LSP-MD：捕捉热振动驱动变构效应的快速计算方法

Molecular Dynamics · 2026-01-16

TS-DAR实用指南：生物分子模拟中的过渡态分析

TS-DAR实用指南：生物分子模拟中的过渡态分析本文信息标题: A Practical Guide to Transition State Analysis in Biomolecular Simulations with TS-DAR（生物分子模拟中的过渡态分析实用指南）作者: Eshani C. Goonetilleke, Bojun Liu, Yue Wu, Michael S. O’Connor, Xuhui Huang 发表时间: 2025年10月31日（接收：2025年8月31日；修订：2025年10月30日；接受：2025年10月31日）单位: Department of Chemistry, Theoretical Chemistry Institute, University of Wisconsin-Madison，美国（美国威斯康星大学麦迪逊分校化学系、理论化学研究所）引用格式: Goonetilleke, E. C., Liu, B., Wu, Y., O’Connor, M. S., & Huang, X. (2025). A Practical Guide to Transition State Analysis in Biomolecular Simulations with TS-DAR. The Journal of Physical Chemistry B. https://doi.org/10.1021/acs.jpcb.5c06097 教程代码: https://github.com/xuhuihuang/ts-dar-tutorials 开源协议: CC-BY 4.0 摘要蛋白质功能所必需的构象变化涉及通过复杂自由能景观中多个短暂、高能态的转变。尽管现有方法如马尔可夫状态模型（MSM）和基于分子动力学（MD）模拟构建的非马尔可夫方法能够有效捕获亚稳态，但它们在识别过渡态方面存在困难。TS-DAR（Transition State Identification via Dispersion and Variational Principle Regularized Neural Networks）是一个计算框架，利用分布外检测（OOD）系统识别特定生物分子构象变化中涉及的所有过渡态。TS-DAR利用深度学习模型将MD模拟中的蛋白构象映射到超球面潜在空间，这种低维表示保留了生物分子构象变化的关键动力学信息。为了区分亚稳态和过渡态，TS-DAR使用VAMP-2和分散损失函数，实现过渡态构象的自动识别。该框架提供了蛋白构象景观的全面视图，促进了药物结合、酶活性和突变效应的研究。核心结论 TS-DAR利用分布外检测技术系统识别生物分子过渡态，解决了传统MSM方法的局限超球面潜在空间映射保留了关键动力学信息，同时实现低维表示 VAMP-2和分散损失函数的组合实现亚稳态与过渡态的自动区分在多个系统中验证，从简单的丙氨酸二肽到复杂的DNA修复蛋白AlkD 在准确性和效率上优于MaxEnt-VAMPNets和MSM-committor等现有方法揭示了蛋白-DNA氢键在AlkD易位速率限制步骤中的关键作用背景蛋白质构象变化是其生物学功能的核心，驱动着酶催化、信号转导和变构调控等关键过程。准确建模这些分子过程对于理解生物分子机制和开发靶向治疗至关重要。马尔可夫状态模型（MSM）和非马尔可夫方法（如quasi-MSM和IGME模型）基于大规模分子动力学模拟构建，是识别亚稳态及表征其转变的强大工具。然而，该领域面临的一个持续挑战是过渡态的识别。过渡态是关键但稀疏填充的构象，定义了分子过程的速率限制步骤。传统的MSM方法虽然能够有效捕获长时间尺度的动力学行为，但由于过渡态在相空间中的瞬态性质和稀疏性，难以直接识别这些高能构象。分布外检测（OOD Detection）的发展为解决这一挑战提供了新思路。OOD检测最初用于提高人工智能在高风险应用（如自动驾驶汽车）中的可靠性，确保系统在遇到不熟悉场景时不会做出错误预测。将OOD检测引入生物分子模拟领域，可以利用其识别异常数据点的能力，将过渡态视为偏离亚稳态分布的分布外构象。关键科学问题本文旨在解决以下核心问题：如何系统识别生物分子过渡态：现有MSM和非马尔可夫方法能有效捕获亚稳态，但对于瞬态、高能的过渡态构象缺乏自动化识别能力如何在低维表示中保留动力学信息：需要一种方法能将高维MD轨迹映射到低维空间，同时保留关键的动力学信息和过渡态特征如何区分亚稳态与过渡态：需要明确的数学框架和损失函数，能够自动区分这两类构象的不同特征如何提供实用的工具和教程：现有方法往往理论性强但缺乏易用的实现，需要提供完整的工作流程和代码教程创新点首次将分布外检测技术应用于生物分子过渡态识别，开创性地将AI安全领域的方法引入计算生物物理超球面潜在空间设计：通过L2归一化约束将特征嵌入映射到固定半径超球面，确保几何一致性 VAMP-2与分散损失的创新组合：VAMP-2损失确保亚稳态内部紧密性，分散损失强制亚稳态间分离端到端自动化框架：从MD轨迹到过渡态识别的完整流程，无需手动调整反应坐标与MSM的无缝集成：TS-DAR状态分配可直接用于构建MSM，提供完整的动力学描述开源教程和代码：提供详细的实现指南和示例代码，降低使用门槛研究内容 TS-DAR框架概述 TS-DAR提供了一个端到端的分析流程，将MD模拟数据转化为对蛋白构象动力学的深入理解。整个工作流程包括四个主要步骤：MD采样、特征化、TS-DAR建模和MSM构建。图1：使用TS-DAR研究蛋白动力学的端到端流程 A.在两个或多个功能构象态之间进行大规模MD模拟 B.选择相关特征来描述研究系统 C.TS-DAR使用神经网络将分子构象映射到超球面潜在空间，提供压缩的低维表示同时保留关键动力学信息 D.TS-DAR的状态分配可用于构建马尔可夫状态模型图2：TS-DAR框架详解 A.TS-DAR使用来自MD轨迹的转移对（$x_t$和$x_{t+\tau}$）作为输入，包含L2归一化层以生成超球面嵌入。Softmax输出用于获得伪状态分配。超球面嵌入和伪状态分配用于估计损失函数。TS-DAR使用包含VAMP-2损失和加权分散损失的组合损失函数优化神经网络 B.L2归一化层将特征嵌入（$\tilde{z}$）限制在半径为$\gamma$的超球面内，得到超球面嵌入（$z$） C.超球面潜在空间。圆圈表示亚稳态数据，星号表示亚稳态中心。实线箭头突出显示状态内紧密性（来自VAMP-2损失），虚线箭头突出显示状态间分散（来自分散损失）工作流程四步骤第一步：MD模拟采样在两个或多个功能构象态之间进行大规模MD模拟，生成覆盖相关构象空间的轨迹数据。模拟需要充分采样亚稳态之间的转变，以捕获过渡态构象。第二步：特征化从高维MD数据中选择能够捕获最相关构象动力学的结构特征。TS-DAR支持两种自动特征选择方法： spectral oASIS：基于变分原理，高效识别最能捕获慢动力学的特征子集 MoSAIC：基于相关性的方法，使用Leiden社区检测算法将相似特征聚类，大簇代表集体运动第三步：TS-DAR建模 TS-DAR的核心是将MD构象映射到结构化潜在空间，实现过渡态的自动检测。第四步：MSM构建使用TS-DAR的状态分配构建马尔可夫状态模型，验证模型的准确性，并预测长时间尺度的动力学行为。 TS-DAR模型架构详解 TS-DAR模型由三个核心部分组成： 1. 编码器神经网络编码器负责处理输入特征，将转移对$(x_t, x_{t+\tau})$从MD轨迹映射到特征嵌入$\tilde{z}$。这里$\tau$是滞后时间，捕获了系统的时间演化信息。 2. L2归一化层 L2归一化约束将所有特征嵌入$\tilde{z}$投影到固定半径$\gamma$的超球面上，得到超球面嵌入$z$： \[z = \gamma \frac{\tilde{z}}{\|\tilde{z}\|}\] 这一步骤至关重要，因为它：确保几何一致性：所有数据点到球心的距离相同便于距离计算：超球面上的距离直接反映构象相似性支持OOD检测：过渡态作为远离亚稳态中心的点更容易识别 3. Softmax输出层 Softmax层输出伪状态分配，用于计算损失函数和估计转移矩阵。损失函数设计 TS-DAR的损失函数结合了两个互补的组分： \[\mathcal{L}_{\text{total}} = \mathcal{L}_{\text{VAMP-2}} + \lambda \mathcal{L}_{\text{dispersion}}\] VAMP-2损失 VAMP-2（Variational Approach for Markov Processes）基于变分原理，最大化潜在空间中慢动力学的捕获能力。该损失函数促进：亚稳态内部的紧密性：同一亚稳态的构象在潜在空间中聚集动力学信息的保留：保持转移矩阵的特征值结构分散损失分散损失强制不同亚稳态在潜在空间中相互分离： \[\mathcal{L}_{\text{dispersion}} = -\sum_{i \neq j} d(c_i, c_j)\] 其中$c_i$和$c_j$是不同亚稳态的中心，$d(\cdot, \cdot)$是超球面上的距离度量。分散损失确保：亚稳态间的清晰边界：不同状态在潜在空间中充分分离过渡态的突出性：位于状态边界的过渡态更容易被识别为OOD点权重参数$\lambda$平衡了两个损失项的贡献，通常需要根据系统特性进行调整。过渡态识别机制 TS-DAR使用OOD分数量化每个构象偏离亚稳态分布的程度： \[\text{OOD}(x) = \min_i d(z(x), c_i)\] 其中$z(x)$是构象$x$的超球面嵌入，$c_i$是第$i$个亚稳态的中心。OOD分数越高，构象越可能是过渡态。通过分析OOD分数的分布，可以设定阈值自动识别过渡态构象。这些构象位于亚稳态之间的边界区域，对应于自由能景观上的鞍点。与MSM的集成 TS-DAR的状态分配可以直接用于构建马尔可夫状态模型：状态定义：TS-DAR自动识别亚稳态和过渡态，提供明确的状态分配转移矩阵估计：基于状态间的转移计数构建转移概率矩阵动力学验证：比较MSM预测的长时间动力学与MD观察数据性质计算：计算平均首次通过时间（MFPT）、平衡分布等动力学性质这种集成方法的优势在于：完整的构象景观描述：同时捕获亚稳态和过渡态自动化流程：无需手动定义反应坐标或committor函数动力学准确性：VAMP-2损失确保慢动力学的正确捕获案例研究：丙氨酸二肽丙氨酸二肽是测试和验证新方法的经典基准系统。该分子具有两个关键二面角（$\phi$和$\psi$），其自由能景观包含多个亚稳态和明确定义的转变路径。图5：丙氨酸二肽的TS-DAR分析结果 TS-DAR应用结果：亚稳态识别：TS-DAR成功识别了$C_{7eq}$、$C_{7ax}$和$\alpha_R$等主要亚稳态过渡态定位：高OOD分数的构象精确定位在自由能景观的鞍点区域动力学一致性：构建的MSM准确预测了状态间的转移速率与已知结果对比：TS-DAR识别的过渡态与基于committor函数的传统方法高度一致案例研究：DNA修复蛋白AlkD AlkD是一种DNA修复蛋白，通过沿双链DNA（dsDNA）易位来识别和修复损伤碱基。理解其易位机制对于认识DNA修复过程至关重要。图6：AlkD在dsDNA上易位的过渡态分析研究发现：过渡态构象特征：TS-DAR识别出易位过程中的多个过渡态构象关键氢键作用：过渡态分析揭示了蛋白-DNA氢键在速率限制步骤中的关键作用机制洞察：氢键的形成和断裂协调了蛋白在DNA上的步进运动与实验一致：识别的过渡态特征与实验观察的易位动力学相符这些新发现是通过传统MSM方法难以获得的，展示了TS-DAR在揭示复杂生物分子机制方面的独特价值。案例研究：绒毛头片蛋白HP35 图7：HP35的4态TS-DAR模型验证曲线 A.VAMP-2损失的验证曲线显示模型成功收敛 B.分散损失验证曲线表明亚稳态中心得到良好分离 HP35是一个35残基的快速折叠蛋白，具有清晰的折叠/去折叠动力学。TS-DAR分析使用了来自D.E. Shaw Research的300微秒全原子轨迹数据。数据处理：特征选择：使用528个C-α原子对之间的距离（残基间隔至少3个）特征降维：通过spectral oASIS从26565个原始特征中筛选出最能捕获慢动力学的特征模型配置：采用4态模型，feat_dim=3，训练30个epoch 关键发现：识别了4个主要的构象态：对应于HP35折叠过程的不同中间态过渡态网络：揭示了状态间复杂的转换网络，包括多条平行路径 MSM验证：Chapman-Kolmogorov检验显示TS-DAR-MSM准确再现长时间尺度动力学隐式时间尺度：ITS图证实模型捕获了系统的主要松弛模式图10：HP35的4态模型及代表性构象 A.超球面潜在空间中的构象分布，虚线指向亚稳态中心 B-E.4个亚稳态的代表性构象叠加（每态5个构象） F.状态2和3之间的过渡态构象 G.状态2和4之间的过渡态构象训练效率：在Apple M3 Mac上，HP35模型训练仅需约20分钟（30个epoch），其中预训练3分钟，完整训练17分钟。案例研究：蛋白磷酸酶2A（PP2A）图14：PP2A的2态TS-DAR模型验证 A.VAMP-2损失验证曲线 B.分散损失验证曲线 PP2A是一种关键的丝氨酸/苏氨酸磷酸酶，参与多种细胞过程。其B56δ调节亚基的突变与智力障碍和癌症相关。TS-DAR揭示了疾病突变如何通过变构途径影响酶活性。研究设计： MD数据：10条100纳秒全原子轨迹，保存间隔10皮秒特征工程：26565个调节亚基B56δ与催化亚基间的成对距离，通过spectral oASIS降维至1000个特征模型配置：2态模型，60个epoch训练（50个预训练+10个正式训练）图15：PP2A活性位点开放度分析 A.活性位点开放距离示意图：C-臂残基572-574质心与$\ce{Mg^{2+}}$离子间距离 B.两个态的活性位点开放度平均概率密度分布，蓝色为State 1（主要采样关闭构象约0.85 nm），红色为State 2（更多采样开放构象）重要发现：微妙构象差异的捕获：TS-DAR成功区分了活性位点的开放和关闭态，尽管这些差异相对微小变构机制洞察：State 1主要采样关闭的活性位点构象（约0.85 nm，对应PDB 8U1X），State 2倾向于更开放的构象激活机制：结果支持通过N-臂和C-臂从全酶核心释放来实现激活的机制疾病突变影响：远离活性位点的疾病相关突变可能通过改变构象集合分布来影响酶活性训练效率：在Apple M3 Mac上约4分钟（60个epoch），预训练2.5分钟，正式训练1.5分钟。实用指南：超参数设置训练TS-DAR模型需要注意以下关键超参数： 1. 随机种子（Random Seed）使用set_random_seed(x)设置随机种子以确保结果可重现。固定种子使得：数据洗牌、权重初始化等随机操作产生相同结果便于问题追踪和性能基准测试不同运行间的结果可比较 2. 验证集比例选择用于验证的数据百分比（通常10-20%）。验证集用于：监控训练过程中的过拟合选择最佳模型检查点评估模型泛化能力 3. 滞后时间（Lag Time）$\tau$ 滞后时间决定了转移对的时间间隔，影响：捕获的动力学时间尺度：较大的$\tau$捕获较慢的过程马尔可夫性：$\tau$应足够大以满足马尔可夫假设数据利用率：$\tau$过大会减少可用的转移对数量建议：从隐式时间尺度分析（ITS）开始，选择慢动力学趋于平稳的$\tau$值。 4. 潜在空间维度潜在空间维度应：足够高以捕获主要的构象自由度足够低以避免过拟合和计算开销通常设置为慢动力学特征值数量的2-3倍 5. 分散损失权重$\lambda$ 权重$\lambda$平衡VAMP-2和分散损失，需要：网格搜索优化：测试不同$\lambda$值（如0.1, 0.5, 1.0, 5.0）监控OOD分数分布：理想的$\lambda$产生明确的亚稳态-过渡态区分系统依赖性：不同系统的最优$\lambda$可能差异较大 6. 批量大小和训练轮次批量大小：平衡计算效率和梯度估计质量（通常1024-8192）训练轮次：监控验证损失，使用早停策略避免过拟合学习率：建议使用学习率衰减策略模型评估与验证图12：HP35的TS-DAR-MSM验证 A.Chapman-Kolmogorov检验比较TS-DAR-MSM预测的状态概率演化（橙色点）与MD轨迹观察值（灰色点）。两者的良好一致性表明TS-DAR-MSM准确捕获了系统的长时间尺度动力学 B.隐式时间尺度图显示三个主要松弛时间尺度随滞后时间的变化验证指标： VAMP-2分数：量化慢动力学捕获能力，分数越高越好 Chapman-Kolmogorov检验：验证MSM的马尔可夫性 OOD分数分布：检查亚稳态和过渡态的清晰分离与已知结果对比：在基准系统上与传统方法比较性能优势：准确性：过渡态识别准确率超过90%（在测试系统上）效率：相比MSM-committor方法，计算时间减少10倍以上鲁棒性：对超参数选择相对不敏感可扩展性：适用于小分子到大型蛋白复合物 Q&A Q1: TS-DAR与传统的committor函数方法有何本质区别？ A1: 自动化程度：TS-DAR无需预先定义反应坐标，而committor方法需要手动选择或优化反应坐标全局视角：TS-DAR一次性识别所有过渡态，committor方法通常只能分析特定转变路径理论基础：TS-DAR基于OOD检测和变分原理，committor基于转移路径采样计算效率：TS-DAR避免了committor计算中的昂贵采样过程，在大型系统上优势明显适用范围：TS-DAR特别适合具有多个过渡态和复杂转变网络的系统 Q2: 超球面嵌入相比普通欧氏空间有什么优势？ A2: 几何一致性：所有数据点到球心距离相同，消除了幅度偏差距离度量意义明确：超球面上的测地距离直接反映构象相似性 OOD检测友好：过渡态作为远离亚稳态中心的点在球面上更容易识别归一化自然性：避免了特征尺度不一致的问题理论保证：变分原理在超球面约束下仍然成立 Q3: 如何为新系统选择合适的超参数？ A3: 滞后时间$\tau$：从ITS分析开始，选择慢过程特征值趋于平稳的$\tau$ 通常从几十到几百皮秒开始尝试验证MSM的马尔可夫性（Chapman-Kolmogorov检验）分散损失权重$\lambda$：先用$\lambda=1.0$训练baseline模型检查OOD分数分布是否能区分亚稳态和过渡态如果区分不明显，增加$\lambda$；如果亚稳态过度分散，减小$\lambda$ 潜在空间维度：从系统主要构象自由度数量的2-3倍开始对于丙氨酸二肽（2个二面角），2-4维即可对于蛋白结构域运动，可能需要10-20维建议策略：使用小规模数据集快速迭代，找到合理范围后在完整数据集上训练 Q4: TS-DAR识别的过渡态如何用于下游分析？ A4: 结构分析：提取高OOD分数的构象进行可视化分析过渡态构象的关键结构特征识别速率限制步骤中的关键相互作用突变效应预测：比较野生型和突变体的过渡态结构分析突变如何改变能垒和过渡态稳定性指导实验设计和理性突变药物设计：识别过渡态特异性结合位点设计稳定或去稳定过渡态的小分子开发过渡态类似物抑制剂动力学建模：构建包含过渡态的详细MSM 计算反应速率和转移路径预测不同条件下的动力学行为 Q5: TS-DAR方法有哪些局限性和适用范围？ A5: 数据需求：需要充分采样过渡态区域的MD轨迹如果过渡态极其罕见，可能需要增强采样方法建议至少观察到几十到几百次转移事件系统大小：原则上可应用于任意大小的系统大型系统需要更多计算资源和训练时间特征选择在大型系统中尤为重要多时间尺度问题：当系统包含多个分离的时间尺度时，单一滞后时间可能不足可能需要多尺度TS-DAR或迭代策略最佳实践：从简单基准系统开始熟悉方法使用多个独立训练验证结果稳定性结合传统方法（如PMF计算）交叉验证关键结论与批判性总结主要贡献 TS-DAR创新性地将OOD检测引入生物分子过渡态识别，提供了系统、自动化的分析框架超球面潜在空间和双损失函数设计实现了亚稳态紧密性与过渡态突出性的平衡在多个基准系统上验证了方法的有效性，从简单模型系统到复杂蛋白揭示了AlkD易位机制中的关键氢键作用，展示了方法在实际研究中的价值提供开源教程和代码，大幅降低了使用门槛，促进方法推广局限性依赖充分采样：TS-DAR需要MD轨迹中包含足够的过渡态构象，对于极其罕见的转变可能需要增强采样超参数敏感性：分散损失权重$\lambda$对不同系统的最优值差异较大，需要一定的调参经验解释性挑战：深度学习模型的黑盒性质使得潜在空间的物理意义不总是直观计算成本：虽然比committor方法更快，但对于非常大的系统，特征选择和模型训练仍需要可观的计算资源验证困难：在缺乏实验或理论基准的新系统上，过渡态识别的准确性难以直接验证未来研究方向基于PDF中讨论的未来展望，TS-DAR可以在多个方向进一步发展： 1. 等变神经网络集成用等变神经网络替代传统特征选择（如spectral oASIS），自动学习分子系统的对称性（旋转、平移不变性）直接处理C-α原子的笛卡尔坐标，消除手动特征工程需求对大型复杂系统尤其有效，能捕获高维非线性特征 2. 增强采样集成 Metadynamics：使用TS-DAR潜在空间的集体变量进行metadynamics，高效采样罕见转变事件 Committor函数构建：将TS-DAR集体变量转换为连接初态和终态的committor函数，提供转变概率估计主动学习：根据OOD分数自适应选择需要额外采样的构象区域 3. 非马尔可夫动力学建模使用TS-DAR状态分配构建广义主方程（GME）模型，更准确描述非马尔可夫动力学捕获记忆效应和非指数松弛行为 4. 实际应用场景抗生素设计：识别细菌RNA聚合酶转录周期中的过渡态，设计特异性抑制剂 PROTAC设计：分析linker-free遭遇复合物的亚稳态界面构象，指导理性linker优化突变效应预测：系统评估疾病相关突变如何改变过渡态能垒和反应路径酶工程：识别酶催化循环的速率限制步骤，通过定点突变降低能垒提升催化效率 5. 技术改进 GPU加速实现以处理超大规模系统开发自动超参数优化工具增强潜在空间的物理可解释性总结：TS-DAR代表了生物分子模拟领域过渡态识别的重大进展。通过将AI安全领域的OOD检测技术引入计算生物物理，TS-DAR提供了一个系统、自动化、高效的框架来识别和分析生物分子功能中的关键过渡态。随着方法的不断完善和应用范围的扩大，TS-DAR有望成为理解蛋白质动力学、药物设计和酶工程的重要工具。

Molecular Dynamics · 2025-11-21

突破时间尺度壁垒：Gen-COMPAS用生成式AI重新定义罕见事件模拟

突破时间尺度壁垒：Gen-COMPAS用生成式AI重新定义罕见事件模拟本文信息标题: Breaking the Timescale Barrier: Generative Discovery of Conformational Free-Energy Landscapes and Transition Pathways 作者: Chenyu Tang, Mayank Prakash Pandey, Cheng Giuseppe Chen, Alberto Megías, François Dehez, Christophe Chipot 发表时间: 2025年10月28日（arXiv预印本）单位: 法国洛林大学、西班牙马德里理工大学、美国芝加哥大学、美国伊利诺伊大学香槟分校引用格式: Tang, C., Pandey, M. P., Chen, C. G., Megías, A., Dehez, F., & Chipot, C. (2025). Breaking the Timescale Barrier: Generative Discovery of Conformational Free-Energy Landscapes and Transition Pathways. arXiv preprint arXiv:2510.24979v1. 摘要蛋白质折叠、变构调节和膜转运等分子转变是生物学的核心，但由于其内在的罕见性，传统分子动力学模拟难以触及。增强采样方法虽然能加速模拟，但成本高昂且依赖于预先定义的任意集合变量，可能引入偏差。本文介绍了Gen-COMPAS框架，一种生成式committor引导的路径采样策略，能够在不预定义变量的情况下、以极低成本重构完整的转变路径。Gen-COMPAS将扩散生成模型（产生物理上合理的中间态）与committor过滤（精准定位过渡态）相结合。通过从这些中间态出发的短时无偏模拟，该方法能在纳秒级时间尺度内快速获得完整的转变路径集合，而传统方法需要数个数量级更多的采样。应用于从微型蛋白到核糖结合蛋白再到线粒体载体的多种体系，Gen-COMPAS高效地获得了committors、过渡态和自由能景观，将机器学习与分子动力学完美统一，为机理研究和实际应用提供了广泛的洞见。核心结论无需预定义集合变量：Gen-COMPAS自动从无偏MD轨迹中学习committor函数，识别控制转变的关键自由度采样效率提升数百倍：Trp-cage蛋白折叠研究中，将所需采样时间从208微秒降至594纳秒，效率提升约350倍同时获得动力学与热力学信息：一次计算即可获得过渡态集合、committor图谱、转变路径和自由能景观适用于复杂异质体系：不同于多数生成模型仅限蛋白质，Gen-COMPAS可显式处理蛋白-配体复合物和膜蛋白体系揭示新的生物学机制：首次明确证明线粒体AAC转运蛋白的闭合态（O-state）是$\ce{ADP^3-}$转运的必经中间态背景分子体系中的罕见转变——如蛋白质折叠、化学反应、分子识别与结合——是理解生物功能的关键。精确绘制过渡态集合（TSE）、主导路径和自由能景观（FEL）对于药物发现、蛋白质设计和催化研究至关重要。传统方法主要依赖两种策略：一是暴力分子动力学（brute-force MD），原则上可提供无偏的动力学视图,但计算资源需求巨大。即使是Anton这类专用超级计算机，其模拟范围仍远短于许多生物学或化学相关过程的真实时间尺度。二是增强采样方法（umbrella sampling、metadynamics、adaptive biasing force等），通过沿预定义的集合变量（CVs）施加偏置来加速采样。然而，这类方法的有效性完全取决于CVs的选择质量——如果关键坐标被遗漏，偏置可能扭曲真实机制或无法有效加速采样。近年来，神经网络的引入带来了数据驱动的CV识别方法（如autoencoder、图神经网络），committor理论的发展也将committor $q$ 确立为最优反应坐标。然而，即使是这些先进方法，仍然严重依赖于增强采样技术。与此同时，生成模型（如Boltzmann generators、MDGen、BioEmu）虽然能直接生成平衡构象，但需要海量训练数据（数百毫秒MD轨迹或大规模实验测量），且生成的构象来自学习的潜空间而非真实分子哈密顿量，可能引入热力学预测的系统性偏差。如何在保持物理严格性的同时实现计算效率，真实捕获罕见分子事件的动力学和热力学，仍是一个艰巨挑战。关键科学问题本文旨在解决以下核心问题：能否在不预定义集合变量的情况下，高效探索罕见事件的转变路径？如何同时获得动力学（committor、过渡态）和热力学（自由能景观）信息？生成模型能否在物理严格性约束下，真正加速罕见事件的采样？该框架能否推广到蛋白-配体复合物、膜蛋白等异质体系？创新点首次将扩散生成模型与committor理论结合：形成强大的迭代反馈循环，自动聚焦于过渡态区域完全无需预定义集合变量：消除CV依赖方法中的系统性偏差原生GPU优化：可扩展至大型生物分子体系适用于异质体系：突破多数生成方法仅限纯蛋白质的局限物理严格性：所有采样源自无偏轨迹，直接操作真实哈密顿量极高效率：采样时间缩短2-3个数量级 Gen-COMPAS框架详解核心工作流程 graph TB subgraph 初始化 direction LR A[亚稳态A和B 1-2 ns无偏MD] A --> B[初始数据集] end subgraph 迭代循环 direction TB C[扩散生成模型 生成中间态] D[Committor预测器 识别q≈0.5的结构] E[目标MD 从A和B收敛至目标] F[短时无偏MD 从过渡态shooting] G[累积数据集] C --> D D --> E E --> F F --> G G --> C end subgraph 下游分析 direction LR H[过渡态识别] I[Committor图谱] J[CCS路径提取] K[自由能景观] end 初始化 --> 迭代循环迭代循环 --> 下游分析 Gen-COMPAS的工作流程可分为三个阶段：阶段1：初始化对两个亚稳态（反应物态A和产物态B）进行极短的无偏MD模拟（1-2 ns）生成初始训练数据集阶段2：迭代优化循环扩散生成模型：基于当前数据集，生成连接A和B的物理合理的中间态构象 Committor预测：在构象空间中学习高维committor函数 $q$，识别 $q \approx 0.5$ 的近过渡态结构（separatrix，即等概率返回A或B的超曲面）目标MD（TMD）：从A和B态出发，收敛至生成的中间态目标 Shooting模拟：从separatrix上的点出发，进行短时无偏MD模拟数据累积：新生成的数据用于下一轮扩散模型和committor预测器的训练阶段3：下游分析识别过渡态集合（TSE）构建投影到任意可解释CVs上的committor图谱提取committor一致性路径（CCS）近似自由能景观（FEL）方法学细节扩散生成模型采用去噪扩散概率模型（DDPM），通过逐步添加噪声将数据分布转化为高斯分布，再通过学习的逆过程生成新样本。关键是该模型在蛋白质构象空间中训练，能够生成物理上合理的中间态。 Committor函数学习 Committor $q(\mathbf{x})$ 定义为从构象 $\mathbf{x}$ 出发、在到达A之前先到达B的概率。$q = 0.5$ 的超曲面（separatrix）精确定义了过渡态集合。Gen-COMPAS直接在笛卡尔空间中用神经网络表示committor，无需预定义CVs。目标MD（TMD）通过在原子坐标上施加调和约束，引导系统从起始态（A或B）收敛至目标构象。这确保生成的中间态能够通过物理上合理的路径到达。 Shooting模拟从separatrix附近的构象出发，进行短时（纳秒级）无偏MD模拟，观察其最终落入A态或B态的概率，以验证和细化committor预测。与现有方法的对比特性传统增强采样纯生成模型 Gen-COMPAS 需要预定义CVs ✓ ✗ ✗ 训练数据需求中等极大（ms级）极小（ns级）物理严格性高（可能有偏）低（潜空间）高（无偏MD）采样效率低-中高（训练后）高动力学信息有限近似精确（committor）异质体系 ✓ ✗（多数） ✓ GPU优化部分 ✓ ✓ 研究内容案例1：Trp-cage快速折叠蛋白研究背景 Trp-cage微型蛋白（20个残基）是快速折叠蛋白研究的经典体系。其简单拓扑和短折叠时间使其成为验证计算方法的理想基准——既小到可以进行全原子模拟，又足够复杂以捕获协同转变和类天然中间态的形成。图2：Gen-COMPAS对Trp-cage快速折叠蛋白的研究图2A：Trp-cage在折叠态、过渡态和解折叠态的代表性结构高亮显示 $\alpha$-螺旋和端到端距离 $d$ 图2B：学习到的committor $q$ 投影到三个集合变量上所有 $\ce{C_\alpha}$ 原子的RMSD $\alpha$-螺旋 $\ce{C_\alpha}$ 原子的RMSD 端到端距离同时显示两条采样路径图2C：自由能景观投影到相同的CVs上左图：Gen-COMPAS结果右图：DESRES微秒级模拟结果（作为对照）图2D：沿路径的归一化CV值和对应的committor值关键发现 Gen-COMPAS成功识别了中间过渡态并重构了折叠自由能景观。估计的自由能差与DESRES的208微秒参考模拟定量一致。Committor和TSE分析揭示了折叠机制呈双路径分叉：路径1：早期螺旋成核，随后核心固化路径2：以中心色氨酸残基周围的疏水塌缩为起始，螺旋形成滞后这一多路径观察与之前的显式溶剂转变路径采样研究一致。效率提升采样时间从208微秒降至594纳秒，效率提升约350倍。更重要的是，这种效率不仅限于平衡自由能景观，还扩展到动力学研究——通过直接学习committor（折叠先于解折叠的概率），Gen-COMPAS无需超长模拟即可表征折叠动力学。案例2：核糖结合蛋白（RBP）的结合伴随折叠研究背景 RBP是一种周质结合蛋白（PBP），对细菌营养摄取至关重要，介导高亲和力核糖识别和递送。与其他PBPs类似，RBP在开放态（apo，无配体）和闭合态（holo，结合配体）之间发生大尺度构象变化。核糖结合与柔性区域的折叠紧密耦合，使RBP成为结合伴随折叠过程的模型体系，展示了局部无序如何驱动分子识别。捕获这一过程在计算上极具挑战性，需要同时采样大尺度运动、局部环结构化和配体对接。图3：Gen-COMPAS对RBP结合伴随折叠过程的研究图3A：RBP-核糖未结合态和结合态的代表性结构、过渡态以及转变路径标注三个集合变量：$d$（配体距离）、hinge（铰链角）、twist（扭转角）图3B：Committor和committor一致性路径揭示两种不同的结合-折叠机制颜色表示committor值图3C：Gen-COMPAS获得的自由能景观左图：$d$ vs hinge 右图：$d$ vs twist 清晰显示结合伴随折叠机制双路径机制 Gen-COMPAS重构了RBP的完整结合伴随折叠路径，捕获了从无序开放态、经部分折叠中间态、到最终核糖结合闭合态的转变，揭示了配体结合与蛋白质折叠的协同相互作用：扭转角重排：两条路径类似，表明运动保守铰链弯曲重排：路径分化路径1：分步诱导契合机制（配体结合先于蛋白闭合）路径2：同步结合与折叠定量热力学与动力学统一 Committor函数直接估计识别了TSE，此时核糖正在与RBP相互作用，但蛋白尚未完全闭合。沿配体位置和域间角度的回溯定义CVs投影的自由能景观，显示出明显的开放和闭合能量盆，由与预期门控运动一致的能垒分隔。这些分析将热力学和动力学统一为结合伴随折叠过程的连贯机理图谱。通过原子水平分辨转变路径并定量其能量学，Gen-COMPAS提供了解剖复杂耦合折叠-结合事件的通用策略，对其他PBPs或本征无序蛋白具有广泛相关性。案例3：线粒体ADP/ATP载体（AAC）的复杂转运机制研究背景 AAC是一种必需的线粒体转运蛋白，通过严格的1:1反转运机制交换胞质$\ce{ADP^3-}$和基质$\ce{ATP^4-}$，维持细胞能量平衡，为生物合成、信号转导和肌肉收缩等过程提供燃料。由于线粒体内膜对核苷酸不可透过，AAC是$\ce{ADP^3-}$的唯一入口，使其功能对氧化磷酸化不可或缺。 AAC通过交替通路机制运作，在胞质开放态（C-state，结合$\ce{ADP^3-}$）和基质开放态（M-state，释放$\ce{ADP^3-}$并结合ATP导出）之间循环。虽然晶体学研究在抑制剂存在下捕获了这些末端态，但并未解释AAC如何在转变过程中防止不受控的核苷酸泄漏。生化和计算研究因此提出了一个瞬态闭合态（O-state），其中$\ce{ADP^3-}$被完全包裹在中央腔内，与膜两侧隔绝。这一状态充当保障机制以确保严格耦合的交换，但直接证据一直难以捕获——常规MD难以捕获这一短寿命中间态，实验数据也大多是间接的。图4：Gen-COMPAS对线粒体ATP/ADP载体（AAC）的研究图4A：holo-AAC（$\ce{ADP^3-}$结合态）的三个亚稳态、过渡态和转变路径从C-state → O-state → M-state的路径 $\ce{ADP^3-}$从O-state到M-state的轴向运动三维CV空间中的构象转变路径和自由能景观能量盆图4B：AAC在膜中的holo-state以及描述构象转变的集合变量 $d_1$：胞质侧门控残基对距离之和 $d_2$：基质侧门控残基对距离之和 $d_3$：$\ce{ADP^3-}$的N6原子到胞质侧门控残基质心的距离图4C：连接三个状态的committor和committor一致性路径图4D：投影到二维的自由能景观（holo-AAC）图4E：apo-AAC的两个亚稳态及其投影到 $d_1$ 和 $d_2$ 的自由能景观明确证实闭合中间态的存在 Gen-COMPAS克服了传统方法的局限，明确展示了$\ce{ADP^3-}$转运过程中闭合中间态的存在。模拟显示转变通过明确定义的 $\text{C} \to \text{O} \to \text{M}$ 路径进行： $\ce{ADP^3-}$首先在胞质开放态中牢固结合通过跨膜螺旋的重排被困在闭合中间态当AAC采用基质开放构象时最终释放到基质中这一路径确认O-state是$\ce{ADP^3-}$导入基质的必经步骤，而非偶然构象。三维集合变量与热力学分析自由能景观投影到三个事后定义的CVs上： $d_1 = \sum (\text{CE-ASP231/CD-LYS32, CZ-ARG137/CE-GLU29, CZ-ARG234/CD-ASP134})$ $d_2 = \sum (\text{CE-LYS95/CD-ASP195, CE-LYS198/CD-ASP291, CE-LYS294/CD-ASP92})$ $d_3 = \text{distance}(\ce{ADP^3-}\ \text{N6}, \text{COM of}\ d_1\ \text{residues})$ 这三个CVs捕获了此转变的关键热力学特征。观察到对应C、O、M态的明显极小值，能垒与AAC预期的螺旋门控运动一致。O-state显现为自由能盆，证明其热力学上足够稳定以充当真正的中间态。 Committor分析的动力学意义 Committor分析进一步确立了O-state的动力学作用：一旦AAC转变到闭合态，向基质开放构象前进的概率占主导，确认这一中间态是向内转运过程的决定性检查点。配体对转运路径的关键作用为进一步探究配体结合与无配体条件下的机理差异，作者对AAC apo-state进行了Gen-COMPAS研究。与holo-state的相对较低能垒（C → O约2.5 kcal/mol，O → M约2 kcal/mol）形成鲜明对比，apo-state中O-state消失，C-M之间出现高达约10 kcal/mol的能垒。这一结果确认apo-state转变在热力学上强烈不利，与之前研究一致，突显了底物在稳定转运兼容路径中的关键作用。生理与疾病意义这些发现不仅阐明了AAC功能的分子机制，还将转运蛋白动力学与线粒体生理和疾病联系起来。由于AAC是生物能量学的核心，$\ce{ADP}$转运功能障碍可能损害ATP生产，导致从线粒体肌病到神经退行性疾病的各种病症。通过确认和定量闭合态，Gen-COMPAS提供了理解突变或抑制剂如何破坏$\ce{ADP^3-}$转运的机理基础，为健康和病理学中线粒体功能的合理调控铺平道路。 Q&A Q1: Gen-COMPAS如何处理大型复杂体系的计算成本？ A1: Gen-COMPAS的计算效率来自三个方面：聚焦采样：通过committor引导，采样集中在过渡态区域（separatrix附近），避免浪费资源探索不相关的构象空间短时模拟：初始化仅需1-2 ns的无偏MD，shooting模拟也在纳秒级，远短于传统方法的微秒-毫秒需求 GPU原生优化：扩散模型和神经网络committor预测器都原生支持GPU加速，可高效处理大型体系（如AAC的膜蛋白体系包含数万原子） Q2: 为什么不能直接使用纯生成模型（如Boltzmann generators）来替代Gen-COMPAS？ A2: 纯生成模型存在以下局限：训练数据需求巨大：通常需要数百毫秒的MD轨迹或大规模实验数据，对罕见事件来说几乎不可能收集物理严格性问题：生成的构象来自学习的潜空间而非真实哈密顿量，可能引入mode imbalance、训练集几何偏好等artifacts 热力学不准确：自由能景观不能与直接MD严格比较，定量预测可靠性有限缺乏动力学信息：无法提供committor等动力学可观测量 Gen-COMPAS通过将生成模型与物理严格的无偏MD结合，克服了这些限制，确保所有结果源自真实分子力学 Q3: Committor函数 $q$ 为什么是最优反应坐标？ A3: 从理论角度，committor $q(\mathbf{x})$ 定义为从构象 $\mathbf{x}$ 出发、在到达反应物态A之前先到达产物态B的概率。它具有以下独特性质：动力学最优性：满足后向Kolmogorov方程，是唯一完全表征转变动力学的单变量定义过渡态：$q = 0.5$ 的超曲面（separatrix）精确定义了过渡态集合，无任意性无需预判：不需要预先知道哪些自由度是慢变量，自动从全维构象空间中提取关键信息路径独立：committor一致性路径（CCS）追踪动力学上有意义的路径，通常不同于最小自由能路径（MFEP） Q4: Gen-COMPAS能否应用于化学反应或材料科学中的相变？ A4: 原则上可以。Gen-COMPAS的核心思想——通过生成模型探索中间态+committor引导识别过渡态——是通用的。只要能定义清晰的两个亚稳态（反应物/产物，或两种相），并能进行短时MD模拟，就可以应用。但需要注意：化学反应可能需要ab initio MD或机器学习力场以准确描述键的断裂/形成相变可能涉及更复杂的多态共存，可能需要扩展到多于两个亚稳态的情况当前实现针对生物分子优化，应用于其他领域可能需要调整生成模型架构 Q5: 如何验证Gen-COMPAS结果的准确性？ A5: 文章采用了多重验证策略：与标准基准对比：Trp-cage结果与DESRES的208微秒全原子模拟定量一致与已知机制一致性：RBP的双路径机制、AAC的O-state中间态均与已有实验和计算证据吻合内部一致性检查：Committor预测通过shooting模拟验证（观察从 $q \approx 0.5$ 构象出发的轨迹确实以约50%概率到达A或B）物理合理性：所有中间态构象通过TMD从真实亚稳态出发可达，确保路径的物理可行性事后CV投影：将结果投影到已知的物理相关CVs上，检查能垒、极小值位置是否合理关键结论与批判性总结主要贡献方法学突破：首次将扩散生成模型与committor理论有机结合，形成完全无CV依赖的罕见事件采样框架效率革命：将罕见事件采样所需时间从微秒-毫秒级降至纳秒级，提升2-3个数量级机理洞见：揭示了Trp-cage的双折叠路径、RBP的双机制结合伴随折叠、AAC的必经闭合中间态等新的生物学机制通用性：适用于蛋白折叠、配体结合、膜转运等多种复杂过程，克服了纯生成模型仅限蛋白质的局限潜在影响药物发现：快速预测配体结合路径和能垒，指导先导化合物优化蛋白质设计：理解折叠路径以设计快速折叠或特定中间态的蛋白转运蛋白研究：解析膜蛋白转运机制，为疾病相关突变和抑制剂设计提供结构基础方法学范式转变：可能改变分子模拟领域对暴力MD和传统增强采样的依赖局限性两态假设：当前框架主要针对两个亚稳态之间的转变，多态体系（如多个折叠中间态、多步骤反应）的扩展尚未充分探索力场依赖：结果准确性仍然受限于底层力场（或机器学习力场）的质量定量精度：作者坦诚指出，生成的自由能景观应作为快速探索工具，对于需要最高定量精度的应用，应以此为起点再用传统增强采样或转变路径采样精修迭代收敛判据：文中未详细讨论如何判断迭代何时收敛，可能需要经验或额外的收敛诊断未来研究方向扩展到多个亚稳态的复杂网络（超越简单两态跃迁）结合更高精度的ab initio MD或神经网络势能面开发自动化的收敛诊断和不确定性量化应用于更大尺度的生物复合物（如病毒衣壳组装、染色质重塑）与实验技术（如单分子FRET、低温电镜）结合验证预测批判性思考正如作者在Discussion中提出的provocative问题：“我们是否正在见证微秒至毫秒级暴力模拟和传统增强采样策略终结的开端？”虽然Gen-COMPAS并非完全取代这些方法，但它确实展示了罕见事件动力学的本质可以在不付出极端计算代价的情况下揭示。然而，我们也需要清醒地认识到，Gen-COMPAS的成功建立在精心选择的基准体系上（Trp-cage、RBP、AAC都是研究充分的体系）。对于全新的、机制完全未知的体系，如何确保生成模型不会产生物理上虽然合理但动力学上不相关的中间态，仍是一个开放性问题。此外，committor的学习本身也需要足够的过渡态附近的采样，如何在初始迭代中高效bootstrap这一过程，可能是方法鲁棒性的关键。尽管如此，Gen-COMPAS无疑代表了分子模拟领域的一次重要进步，将机器学习的生成能力与物理模拟的严格性完美结合，为我们理解生命分子的动态过程打开了新的窗口。

Molecular Dynamics · 2025-11-17

靶向分子动力学（TMD）：用RMSD约束引导蛋白质构象转变

Molecular Dynamics · 2025-10-11

自适应偏置力（ABF）方法详解

自适应偏置力（ABF）方法详解一、ABF方法的基本原理自适应偏置力（Adaptive Biasing Force, ABF）是一种用于计算自由能曲面（PMF）的增强采样方法。它的核心思想是：通过实时计算并施加一个抵消系统平均力的偏置力，使分子能够在反应坐标上自由扩散，从而加速采样。基本方程对于一个集合变量（collective variable, CV）$\xi$，系统在 $\xi$ 方向上受到的瞬时力为 $F(\xi)$。ABF方法通过累积统计，估算出在 $\xi$ 处的平均力 $\langle F(\xi) \rangle$： \[\langle F(\xi) \rangle = -\frac{\mathrm{d}A(\xi)}{\mathrm{d}\xi}\] 其中 $A(\xi)$ 是沿着 $\xi$ 的自由能（PMF）。 ABF的策略：在模拟过程中，实时施加一个偏置力 $F_{bias}(\xi) = -\langle F(\xi) \rangle$，使得分子在 $\xi$ 方向上受到的净力接近零，从而能够自由地在整个 $\xi$ 范围内扩散。瞬时力的计算：从原子力到集合变量的投影关键问题：MD引擎（如NAMD、GROMACS）计算的是原子间的相互作用力 $\mathbf{F}_i$（作用在每个原子 $i$ 上），但ABF需要的是沿着集合变量 $\xi$ 的广义力 $F(\xi)$。如何将原子力转换为CV方向的力？答案：通过链式法则投影。集合变量 $\xi$ 通常是原子坐标 ${\mathbf{r}_i}$ 的函数，即 $\xi = \xi(\mathbf{r}_1, \mathbf{r}_2, \ldots, \mathbf{r}_N)$。瞬时力通过以下公式计算： \[F(\xi) = -\sum_{i=1}^{N} \mathbf{F}_i \cdot \frac{\partial \xi}{\partial \mathbf{r}_i}\] 物理意义： $\frac{\partial \xi}{\partial \mathbf{r}_i}$ 是CV对第 $i$ 个原子坐标的梯度，表示该原子沿哪个方向运动会增加 $\xi$ 的值 $\mathbf{F}_i \cdot \frac{\partial \xi}{\partial \mathbf{r}_i}$ 是原子 $i$ 受到的力在CV方向上的投影分量负号是因为力的定义（$\mathbf{F} = -\nabla U$）具体例子：在本文中，CV是小分子沿膜法线（z轴）的位置，即 $\xi = z_{molecule}$。此时： $\frac{\partial \xi}{\partial \mathbf{r}_i} = (0, 0, 1)$ 只有z分量非零 $F(\xi) = -F_{i,z}$ 只需提取分子受力的z分量实际实现：每个MD时间步，MD引擎计算所有原子受到的力 ${\mathbf{F}_i}$ Colvars模块（NAMD）或相应的插件（GROMACS）实时计算：当前的CV值 $\xi(t)$ CV的梯度 ${\partial\xi/\partial\mathbf{r}_i}$ 瞬时广义力 $F(\xi,t)$ 累积到直方图：将 $F(\xi,t)$ 加到对应 $\xi$ 网格点的累积和中计算平均力：$\langle F(\xi) \rangle = \frac{1}{N_{samples}(\xi)} \sum_{t:\xi(t)\approx\xi} F(\xi,t)$ 施加偏置：在下一个时间步，对相关原子施加偏置力 $\mathbf{F}_{bias,i} = -\langle F(\xi) \rangle \cdot \frac{\partial \xi}{\partial \mathbf{r}_i}$ 技术细节： ABF使用分层网格将CV空间离散化（如每0.01 nm一个网格点）为避免初期统计不准确，通常设置最小采样阈值（如每个网格点至少100次访问）才开始施加偏置力偏置力的施加使用渐进式缩放（ramp），从0逐渐增加到1，避免非平衡效应自由能的恢复模拟结束后，通过对累积的平均力进行积分，即可恢复自由能曲面： \[A(\xi) = A(\xi_0) - \int_{\xi_0}^{\xi} \langle F(\xi') \rangle \mathrm{d}\xi'\] 二、ABF的窗口策略与边界处理为什么需要分窗口？虽然理论上ABF可以在整个反应坐标范围内一次性进行（全局ABF），但在实际应用中，当自由能曲面存在高能垒时，全局ABF会遇到严重的采样问题：能垒区域采样不足：分子很难跨越高能垒区域，导致这些区域的平均力估计不准确收敛极慢：即使施加了偏置力，分子在能垒区域的停留时间仍然很短，需要极长的模拟时间才能充分采样解决方案：将整个反应坐标范围划分为多个重叠的窗口（stratification），在每个窗口内独立进行ABF采样，最后将各窗口的PMF拼接起来。窗口的定义每个窗口由以下参数定义：窗口范围 $[\xi_{min}, \xi_{max}]$：CV允许的取值范围窗口宽度：$\Delta\xi = \xi_{max} - \xi_{min}$（本文中为0.4 nm）窗口中心：$\xi_{center} = (\xi_{min} + \xi_{max})/2$ 相邻窗口的间隔：中心点之间的距离（本文中为0.1 nm）例如，在本文中：窗口1：$[-0.2, +0.2]$ nm，中心在 0 nm 窗口2：$[-0.1, +0.3]$ nm，中心在 +0.1 nm 窗口3：$[0.0, +0.4]$ nm，中心在 +0.2 nm … 边界的处理方式 ABF方法对窗口边界的处理与umbrella sampling有本质区别： 1. 无强制约束的边界 ABF不在窗口边界施加强制约束势。当CV的值 $\xi$ 处于窗口范围 $[\xi_{min}, \xi_{max}]$ 内时：正常施加偏置力：$F_{bias}(\xi) = -\langle F(\xi) \rangle$ 正常采样和累积统计：该位置的构象被记录用于平均力的估算当 $\xi$ 超出窗口范围时：停止施加偏置力：不再对系统施加ABF偏置停止采样：该位置的构象不被记录模拟继续运行：系统仍然正常演化，只是不参与当前窗口的统计 2. 可选的软约束势（wall potential）为了防止分子过度偏离窗口范围，可以在边界外侧添加一个软约束势（也称为wall potential或restraining potential）： \[U_{wall}(\xi) = \begin{cases} \frac{k}{2}(\xi - \xi_{max})^2 & \text{if } \xi > \xi_{max} + \delta \\ 0 & \text{if } \xi_{min} - \delta \leq \xi \leq \xi_{max} + \delta \\ \frac{k}{2}(\xi - \xi_{min})^2 & \text{if } \xi < \xi_{min} - \delta \end{cases}\] 其中： $k$ 是弹簧常数（通常为10-100 kcal/mol/Å²） $\delta$ 是缓冲区宽度（通常至少为一个网格间距）关键特点：约束势的作用范围应比窗口范围更宽（$\delta > 0$），确保在窗口边界处没有突变约束势是柔和的（软约束），不会强制将分子”锁死”在某个位置与Umbrella Sampling的对比特性 ABF Umbrella Sampling 窗口定义定义边界范围 $[\xi_{min}, \xi_{max}]$ 定义中心点 $\xi_0$ 约束方式无强制约束（或软约束）强制谐振子势 $\frac{k}{2}(\xi-\xi_0)^2$ 分子运动在整个窗口内自由扩散被”拴”在中心点附近，受弹簧限制偏置力动态调整，实时抵消平均力静态谐振子势后处理不需要，直接积分平均力得PMF 需要WHAM等方法去除偏置先验知识不需要知道自由能形状需要预估PMF形状来设置弹簧常数窗口重叠不强制要求（但推荐）必须重叠，否则WHAM无法拼接三、窗口的拼接与PMF的构建重叠区域的作用虽然ABF在理论上不强制要求窗口重叠（因为平均力是连续的），但在实践中高度推荐使用重叠窗口，原因如下：提高统计精度：重叠区域被两个窗口同时采样，提供了交叉验证平滑过渡：减少拼接时的不连续性检测采样质量：如果两个窗口在重叠区域的PMF差异很大，说明采样不充分拼接算法详解 ABF窗口拼接的核心挑战在于：每个窗口独立模拟得到的PMF只是相对值（积分常数未定），需要通过重叠区域将它们”对齐”到同一个能量基准上。步骤1：对每个窗口内的平均力进行积分对于第 $i$ 个窗口（范围 $[\xi_i^{min}, \xi_i^{max}]$），从下边界开始积分平均力： \[A_i(\xi) = -\int_{\xi_i^{min}}^{\xi} \langle F_i(\xi') \rangle \mathrm{d}\xi', \quad \xi \in [\xi_i^{min}, \xi_i^{max}]\] 注意：这里人为设定 $A_i(\xi_i^{min}) = 0$，所以 $A_i(\xi)$ 只是窗口内的相对PMF 积分通常使用数值方法（如梯形法则或辛普森法则）如果平均力在某些点采样不足，可能需要平滑处理（如样条插值）步骤2：在重叠区域对齐相邻窗口对于相邻的窗口 $i$ 和 $i+1$，它们的重叠区域是 $[\xi_{i+1}^{min}, \xi_i^{max}]$。在这个区域内，两个窗口都提供了PMF估计：$A_i(\xi)$ 和 $A_{i+1}(\xi)$。目标：找到一个偏移常数 $\Delta A_i$，使得 $A_i(\xi) + \Delta A_i \approx A_{i+1}(\xi)$ 在重叠区域内尽可能一致。方法1：简单平均法 $\Delta A_i = \frac{1}{N_{overlap}} \sum_{\xi \in overlap} [A_{i+1}(\xi) - A_i(\xi)]$ 方法2：加权最小二乘法（推荐）考虑到不同位置的采样质量不同，使用加权最小二乘： \[\Delta A_i = \arg\min_{\Delta} \sum_{\xi \in overlap} w(\xi) [A_{i+1}(\xi) - A_i(\xi) - \Delta]^2\] 其中权重 $w(\xi)$ 通常取为该点的采样次数：$w(\xi) = \min(N_i(\xi), N_{i+1}(\xi))$，确保采样好的区域有更高的权重。方法3：基于平均力的直接拼接更精确的方法是直接在重叠区域比较平均力，而非PMF： \[\Delta A_i = -\int_{\xi_{i+1}^{min}}^{\xi_i^{max}} [\langle F_{i+1}(\xi') \rangle - \langle F_i(\xi') \rangle] \mathrm{d}\xi'\] 这种方法对噪声更鲁棒，因为它利用了原始的平均力数据。步骤3：全局拼接从第一个窗口开始，逐步累积偏移量，构建全局PMF： \[A(\xi) = \begin{cases} A_1(\xi) & \text{if } \xi \in [\xi_1^{min}, \xi_1^{max}] \\ A_2(\xi) + \Delta A_1 & \text{if } \xi \in [\xi_2^{min}, \xi_2^{max}] \\ A_3(\xi) + \Delta A_1 + \Delta A_2 & \text{if } \xi \in [\xi_3^{min}, \xi_3^{max}] \\ \vdots \\ A_i(\xi) + \sum_{j=1}^{i-1} \Delta A_j & \text{if } \xi \in [\xi_i^{min}, \xi_i^{max}] \end{cases}\] 在重叠区域的处理：对于重叠区域 $[\xi_{i+1}^{min}, \xi_i^{max}]$，可以：选择其一：只使用窗口 $i$ 或窗口 $i+1$ 的数据加权平均（推荐）： $A(\xi) = \frac{w_i(\xi) \cdot [A_i(\xi) + \sum_{j=1}^{i-1}\Delta A_j] + w_{i+1}(\xi) \cdot [A_{i+1}(\xi) + \sum_{j=1}^{i}\Delta A_j]}{w_i(\xi) + w_{i+1}(\xi)}$ 其中 $w_i(\xi) = N_i(\xi)$ 是窗口 $i$ 在 $\xi$ 处的采样次数步骤4：质量检查拼接完成后，应检查：连续性：相邻窗口的PMF在重叠区域是否平滑连接一致性：重叠区域内两个窗口的PMF差异是否小于统计误差（通常 < 0.5 kcal/mol）平均力一致性：重叠区域内 $\langle F_i(\xi) \rangle$ 和 $\langle F_{i+1}(\xi) \rangle$ 是否接近与WHAM的对比： ABF拼接：简单、直接，只需在重叠区域对齐PMF，不需要迭代求解 WHAM：用于umbrella sampling，需要迭代求解自洽方程，计算复杂度更高，但在窗口重叠较少时更稳定四、ABF的优势与局限优势无需先验知识：不需要预先知道自由能曲面的形状高效采样：在能垒高的区域，ABF比umbrella sampling更高效无后处理：不需要WHAM等复杂的后处理方法局限初期采样问题：在模拟初期，平均力估计不准确，需要设置一个最小采样阈值（如每个网格点至少100次访问）才开始施加偏置隐藏能垒：如果正交于CV的自由度存在高能垒，ABF可能采样不充分几何约束的影响：当CV与几何约束或其他CV耦合时，需要使用扩展ABF（extended ABF, eABF）来正确处理五、主流MD软件中的ABF实现 5.1 NAMD中的ABF 实现方式：ABF在NAMD中通过Colvars模块（Collective Variables Module）实现，是NAMD内置的官方支持方法。基本使用流程：定义集合变量：在配置文件中定义CV（如距离、角度、二面角、RMSD等） colvar { name myDistance distance { group1 { atomNumbers 1 2 3 } group2 { atomNumbers 10 11 12 } } } 启用ABF：配置ABF参数 abf { colvars myDistance fullSamples 200 # 开始施加偏置前的最小采样数 historyfreq 50000 # 输出频率 writeTISamples yes # 输出统计数据 } 运行模拟：NAMD自动计算瞬时力、累积平均力并施加偏置支持的集合变量类型： distance：原子间距离 angle、dihedral：键角和二面角 rmsd：相对参考结构的RMSD gyration：回旋半径 eigenvector：沿主成分的投影输出文件： .pmf：PMF曲线数据 .count：每个网格点的采样次数 .grad：平均力数据参考资源： NAMD官方ABF教程：https://www.ks.uiuc.edu/Training/Tutorials/namd/ABF/ Colvars参考手册：https://colvars.github.io/colvars-refman-namd/ 5.2 GROMACS中的ABF 实现方式：GROMACS本身不直接支持ABF，但有以下几种替代方案：方案1：GROMACS + PLUMED（不推荐用于ABF） PLUMED是一个通用的增强采样插件，支持多种MD引擎局限：PLUMED不计算二阶导数，只能实现基于一阶导数的简化ABF版本 ABF并非PLUMED的原生方法，需要自行用C/C++实现方案2：GROMACS + SSAGES（推荐用于ABF） SSAGES（Software Suite for Advanced General Ensemble Simulations）提供了完整的ABF实现使用流程：使用GROMACS工具准备输入文件（拓扑、坐标）编写SSAGES的JSON配置文件定义CV和ABF参数使用gmx_ssages或gmx_mpi运行模拟文档：https://ssagesproject.github.io/ 方案3：GROMACS原生AWH方法（推荐替代） AWH（Accelerated Weight Histogram）是GROMACS 2018及以后版本的原生自适应偏置方法原理类似ABF：通过自适应调整偏置势来加速采样并计算PMF 优势： GROMACS原生支持，无需外部插件性能优化好，与GROMACS集成度高文档完善基本使用： pull = yes pull-ncoords = 1 pull-coord1-type = umbrella pull-coord1-geometry = distance pull-coord1-groups = 1 2 awh = yes awh-nstout = 1000 awh-nbias = 1 awh1-ndim = 1 awh1-dim1-coord-index = 1 参考文档：https://manual.gromacs.org/current/reference-manual/special/awh.html 推荐方案对比：方案优势劣势适用场景 SSAGES 完整ABF实现需要额外编译安装需要严格使用ABF算法 AWH 原生支持、性能好与标准ABF略有差异大多数自适应偏置应用 PLUMED 通用性强、功能多 ABF支持有限使用其他增强采样方法 5.3 其他MD软件 LAMMPS：通过Colvars模块支持ABF（与NAMD共用） Amber：通过PLUMED插件支持有限的ABF功能 OpenMM：通过Colvars或PLUMED插件支持总体建议：如需使用标准ABF方法，NAMD是首选（原生支持，文档完善） GROMACS用户建议使用AWH方法（原生、高效）或SSAGES（标准ABF）对于多维复杂CV或需要与其他增强采样方法结合，考虑使用PLUMED

Molecular Dynamics · 2025-10-09

分子动力学揭示药物靶点变构通信路径：从动态网络到功能调控

title: “MDPath：追踪蛋白质中的“悄悄话”——用分子动力学揭示药物靶点（如GPCRs）的变构通信路径” date: “2025-10-02” tags: [molecular-dynamics, sampling-and-analysis] — MDPath：追踪蛋白质中的“悄悄话”——用分子动力学揭示药物靶点（如GPCRs）的变构通信路径本文信息标题: MDPath：通过分子动力学模拟揭示药物靶点的变构通讯路径作者: Niklas Piet Doering, Marvin Taterra, Marcel Bermúdez, and Gerhard Wolber 发表时间: 2025年9月23日 (Accepted) 单位: 柏林自由大学生物、化学与药学系 (德国)，明斯特大学药物与医药化学研究所 (德国) 引用格式: Doering, N. P., Taterra, M., Bermúdez, M., & Wolber, G. MDPath: Unraveling Allosteric Communication Paths of Drug Targets through Molecular Dynamics Simulations. Journal of Chemical Information and Modeling. Published online September 23, 2025. https://doi.org/10.1021/acs.jcim.5c01590 源代码: https://github.com/wolberlab/mdpath 摘要理解蛋白质中的变构通讯对于基于结构的理性药物设计仍然是一个关键挑战。我们在此推出MDPath，一个用于分析分子动力学模拟中变构通讯路径的Python工具包，其核心是基于归一化互信息（NMI）的分析。我们以β₂-肾上腺素能受体、腺苷A₂A受体和μ-阿片受体为模型系统，展示了MDPath识别已知及新型GPCR变构机制的能力。该工具包揭示了β₂-肾上腺素能受体和MOR中配体特异性的变构效应，阐明了蛋白质-配体相互作用如何驱动构象变化。通过对ABL1激酶与变构和正构抑制剂复合物的分析，证明了该方法的广泛适用性。最终，MDPath为绘制蛋白质内部的变构通讯提供了一个开源框架，推动了基于结构的药物设计。背景变构（Allostery）是生物学中最基本的调控原则之一，它描述了一种“隔山打牛”的现象：蛋白质上一个位点的扰动（如配体结合或氨基酸突变）能够引起远处另一个功能位点的活性发生改变。这种远程调控使得药物分子不必直接作用于蛋白质的活性中心，而是可以通过结合在一个全新的“变构口袋”，来精细地调节蛋白质的功能，这为开发高选择性、低副作用的药物提供了巨大机遇。GPCRs、激酶等许多重要药物靶点都受到变构调控。然而，识别连接这两个远距离位点的“通讯线路”是一个巨大的挑战。这些线路并非静态的物理连接，而是由蛋白质内部残基间动态的、协同的运动所构成的复杂网络。静态的晶体结构往往无法揭示这些隐藏的动态信息，因此，分子动力学（MD）模拟成为捕捉蛋白质动态行为、研究变构机制不可或缺的工具。近年来，虽然涌现出多种用于分析MD轨迹以识别变构网络的计算工具，但它们大多关注于蛋白质整体的通讯网络，难以精确地分离出由特定配体结合所诱导的信号通路。此外，许多工具的设置复杂或并非开源，限制了其在药物研发领域的广泛应用。因此，亟需一个易于使用、开源且能系统性地、定量地描绘配体特异性变构路径的工具。关键科学问题如何从分子动力学模拟的海量数据中，系统性地、自动化地识别并可视化连接药物结合位点与功能远端位点的变构通讯路径？我们能否开发一个通用工具，不仅能确认已知的变构机制（如GPCR中的保守“微开关”），还能揭示配体特异性的调控网络（如激动剂和拮抗剂引发的不同信号通路），并为实验中观察到的突变效应提供合理的动力学解释？创新点发布MDPath开源工具包：提供了一个完整的、从MD轨迹分析到三维可视化的Python工具包，用于系统性地研究蛋白质变构通讯，其代码已在GitHub上开源。基于归一化互信息（NMI）的路径识别：采用NMI来量化残基间动态运动的相关性，并结合图论算法（Dijkstra）来寻找“信息流”最优的路径，为变构分析提供了数学上严谨且物理上直观的方法。配体特异性路径分析：实现了从特定配体接触残基出发追踪通讯路径的功能，能够清晰地区分不同配体（如激动剂与拮抗剂）引发的不同变构信号网络。广泛的验证与应用：在GPCRs和激酶这两大类重要药物靶点上成功验证了该方法，不仅重现了已知的保守变构基序，还为实验突变数据提供了新的机理见解。研究内容分子动力学模拟方法体系构建与参数化：研究使用了多个GPCR体系和ABL1激酶体系。GPCR结构来源于PDB数据库，包括激动剂结合态（β2：7DHI，A2A：2YDO，MOR：8EFQ）和拮抗剂/反向激动剂结合态（β2：5JQH，A2A：5MZP，MOR：7UL4），ABL1激酶结构为8SSN。所有体系使用MOE 2022.2进行预处理，包括缺失环区建模、突变回归野生型序列、添加缺失原子等。模拟软件与力场： GPCR体系：使用OpenMMDL进行体系构建，OpenMM进行MD模拟 ABL1体系：使用CHARMM GUI进行体系构建力场选择：蛋白质使用AMBER14SB力场，脂质使用Lipid21力场，配体使用GAFF2力场（ABL1体系中阿西米尼使用OpenFF）溶剂模型：TIP3P水模型，0.15 M NaCl离子浓度模拟参数：所有体系均进行能量最小化和0.5 ns平衡后，在NPT系综下运行3个独立的200 ns生产模拟。温度控制在300 K（Langevin动力学），压强控制在1.0 atm，时间步长2 fs，每个重复记录1000帧轨迹用于后续分析。核心方法论深度解析：MDPath的工作原理与流程 MDPath的核心思想是将蛋白质看作一个信息传递网络，利用MD模拟捕捉其动态行为，再通过信息论和图论的工具来寻找信息传递效率最高的“高速公路”。图5：MDPath用于变构通讯路径检测的主要工作流程。工作流程分为三个主要阶段：输入阶段接收MD模拟轨迹文件（PDB拓扑和DCD轨迹），可选择性添加配体相互作用位点等参数；分析阶段首先计算残基主链二面角运动，然后计算归一化互信息矩阵，接着构建网络图并使用Dijkstra算法寻找最大NMI路径，最后进行层次聚类识别核心通路；可视化阶段生成多种格式的输出文件，包括NGL view（Jupyter笔记本）、PyMOL脚本和STL文件（用于ChimeraX等软件）。 1. 相关性分析：从原子运动到信息网络为什么选择二面角？ MDPath选择监测每个残基的主链二面角（$\phi, \psi$）的动态变化，而不是Cα原子的笛卡尔坐标。这是一个关键的方法学选择。因为笛卡尔坐标会受到蛋白质在模拟盒子中整体平动和转动的影响，直接计算坐标相关性会引入大量虚假的、无物理意义的噪声。而二面角是内坐标，它只描述了肽链局部的扭转运动，与分子的整体运动无关。因此，基于二面角计算出的相关性更能反映蛋白质内部真实的构象变化和信息传递，信噪比更高。如何量化“通讯”强度？——互信息与NMI MDPath采用信息论中的归一化互信息（Normalized Mutual Information, NMI）来量化任意两个残基（X和Y）之间的“通讯强度”。首先，计算两个残基二面角运动之间的互信息（Mutual Information, MI）： $MI(X,Y)=\sum_{x}\sum_{y}P(x,y)\log_{2}\left(\frac{P(x,y)}{P(x)\cdot P(y)}\right)$ 公式的通俗解释互信息衡量了知道一个变量后，另一个变量不确定性减少的程度，可以理解为两个变量之间非线性相关性的量度。 $P(x,y)$ 是联合概率分布，表示残基X处于状态x（某个二面角角度范围）且同时残基Y处于状态y的概率。 $P(x)$ 和 $P(y)$ 是边缘概率分布，分别表示X处于状态x和Y处于状态y的概率。如果X和Y的运动完全独立，那么 $P(x,y) = P(x) \cdot P(y)$，比值为1，$\log_2(1)=0$，MI为0。如果X和Y的运动高度相关，那么 $P(x,y)$ 会远大于 $P(x) \cdot P(y)$，比值大于1，$\log_2$项为正，MI值就高。然后，使用每个残基自身的熵（Entropy） $H(X)=-\sum P(x)\log_{2}(P(x))$ 对MI进行归一化，得到NMI： $NMI(X, Y) = \frac{MI(X, Y)}{\sqrt{H(X)\cdot H(Y)}}$ 公式的通俗解释 NMI通过除以两个残基各自信息熵的几何平均值，消除了变量自身复杂性的影响。这使得NMI的取值范围被限定在0（完全无关）到1（完全相关）之间。一个高的NMI值意味着两个残基在动态运动上是高度协同的，即使它们在空间上相距很远，也表明它们之间存在一条有效的“通讯”通路。 2. 基于图的路径分析：寻找最优通讯路径计算出所有残基两两之间的NMI值后，MDPath将蛋白质抽象成一个网络图（Graph）。它将每个氨基酸残基视为一个节点（node），并在空间上邻近（< 5 Å）的残基之间创建边（edge）。关键的一步是如何利用Dijkstra算法。Dijkstra算法是图论中一个经典的最短路径算法，它寻找的是图中两点之间权重之和最小的路径。然而，我们的目标是寻找累积NMI值最大（即信息流最强）的路径。为了利用Dijkstra算法，MDPath进行了一个巧妙的转换：它将每条边的权重（weight）定义为与NMI值成反比的量（例如 $w = 1 - NMI$）。这样，NMI值越高（通讯越强），边的权重就越小。因此，在这个权重被“反转”的图中寻找“最短路径”，就等价于在原始概念中寻找“信息量最大的路径”。通过对所有可能的残基对运行该算法并筛选，MDPath便可描绘出蛋白质内部主要的变构通讯网络。 graph TD subgraph "输入阶段" direction LR A1["**MD模拟轨迹** 拓扑文件PDB"] A2["**轨迹文件** DCD格式"] A3["**可选参数** 配体相互作用位点 分析参数设置"] end subgraph "相关性分析阶段" direction LR B1["计算所有残基 主链二面角φψ轨迹"] B2["计算残基对间 归一化互信息NMI矩阵"] B1 --> B2 end subgraph "路径分析阶段" direction LR C1["构建网络图 残基为节点NMI为边权重"] C2["Dijkstra算法 寻找最大NMI路径"] C3["层次聚类 识别核心通路"] C1 --> C2 --> C3 end subgraph "可视化输出阶段" direction LR D1["**NGL view** Jupyter交互式"] D2["**PyMOL脚本** 结构渲染"] D3["**STL文件** ChimeraX等软件"] end A1 --> B1 A2 --> B1 A3 --> B1 B2 --> C1 C3 --> D1 C3 --> D2 C3 --> D3 结果与分析 1. 模拟体系的质量控制：确保动力学轨迹的可靠性图S3-S5：激动剂结合的GPCR在200 ns模拟过程中的A100激活指数变化。 A100激活指数的计算原理：A100是一个专为A类GPCR设计的通用激活指数，基于五个关键的跨膜螺旋间距离计算得出。该指数通过机器学习方法训练，使用了大量微秒级分子动力学模拟数据和268个已发表的X射线晶体结构进行验证。A100指数的分类准确性在二态模型中达到94%（活性态）和99%（非活性态），在三态模型（包括中间态）中对活性态、中间态和非活性态的准确性分别为63%、81%和89%。在分析通讯路径之前，必须确保MD模拟本身是可靠的，即蛋白质在模拟过程中保持在预期的功能状态（活性态或非活性态）。作者使用A100激活指数来监测GPCR的构象状态（分数 > 0表示活性态，分数 < 0表示非活性态）。补充材料中的图S3-S5显示，在所有激动剂结合的体系中，A100分数在200 ns的模拟时长内基本都保持在0以上，表明模拟轨迹很好地维持了受体的活性构象，为后续的路径分析提供了可靠的数据基础。 2. 验证：识别GPCR中的保守变构“微开关” 图1：(A) 沙丁胺醇结合的活性态β₂-肾上腺素能受体的完整路径图。(B) 卡拉洛尔结合的非活性态β₂-肾上腺素能受体的完整路径图。(C) 热图显示了在所有三个模拟重复的前500条路径中，A类GPCR保守基序残基的参与情况。图中蓝色和紫色路径表示变构通讯路径，路径的粗细反映通讯强度。子图(D-H)详细展示了特定基序的路径：蓝色路径穿过CWxP基序(D)和PIF基序(E)，橙色残基标记关键基序位点。在非活性态中，蓝色路径通过NPxxY基序(F)和DRY基序的离子锁结构(G,H)。热图计算方法：图1C的热图统计了前500条最强通讯路径中每个保守基序残基的出现次数。对于每个基序（如CWxP、PIF、NPxxY、DRY），计算该基序内所有残基在路径中的参与频率，然后取该基序内任一残基的最大出现频率作为该基序的代表值。这种计算方式能够量化不同功能状态下各个保守”微开关”基序在变构通讯网络中的重要性。热图使用对数标度以更清晰地显示频率差异，颜色越深表示该基序在相应条件下的参与度越高。 GPCR的激活过程依赖于几个保守的氨基酸基序（”微开关”）的协同运动。MDPath的分析结果与已知的生物学机制高度吻合。在活性态受体(A)中，可以看到从细胞外域延伸到细胞内域的蓝色路径。非活性态受体(B)显示不同的路径模式。如图1C热图所示，在激动剂结合的活性态受体中，与激活相关的CWxP和PIF基序在通讯路径中的出现频率非常高。相反，在反向激动剂结合的非活性态受体中，与稳定非活性态相关的NPxxY和DRY基序则占据了主导地位。 3. 解释：为实验突变数据提供机理支撑图2：(A) 腺苷结合的腺苷A₂A受体中，从T88到W246的路径。(B) DAMGO结合的μ-阿片受体中，通过关键枢纽Y328的路径。图中蓝色路径表示变构通讯路径，橙色残基标记关键位点，黄色分子为配体。在A₂A受体(A)中，蓝色路径连接T88³·³⁶（橙色）到激活开关W246⁶·⁴⁸（橙色），展示从TM3到CWxP基序的直接变构通讯，解释了T88突变导致受体活性降低的机理。在μ-阿片受体(B)中，蓝色路径汇聚于关键枢纽残基Y328⁷·⁴³（橙色），该残基位于NPxxY基序上方，作为路径分布中心控制向细胞内结构域的信号传递。实验表明，在A₂A受体中将T88突变会显著降低受体活性。MDPath的分析（图2A）首次发现了一条从T88直达激活开关CWxP基序的变构路径，为该实验现象提供了清晰的机理解释。同样，对于μ-阿片受体（MOR），MDPath也发现Y328是一个关键的路径“枢纽”（hub）（图2B），与其实验功能的重要性相符。 4. 洞察：绘制配体特异性的通讯网络图3：β₂-肾上腺素能受体中的配体特异性路径。(A) 激动剂沙丁胺醇结合的活性态中的路径集群。(B) 反向激动剂卡拉洛尔结合的非活性态中的路径集群。图中展示了两种不同的变构路径集群：蓝色和红色路径代表两个主要的通讯集群，路径粗细反映通讯强度。黄色分子为配体（沙丁胺醇或卡拉洛尔），橙色残基标记参与路径的关键位点。在激动剂沙丁胺醇结合的活性态(A)中，路径主要汇聚到激活相关的PIF基序，显示出典型的激活信号传递模式。在反向激动剂卡拉洛尔结合的非活性态(B)中，路径模式完全不同，主要连接到稳定非活性态的NPxxY基序。值得注意的是，N312⁷·³⁹在两种状态下都不是主要路径的组成部分，表明其主要作用可能是配体结合而非功能调控。 5. 方法的稳健性与拓展应用模型完整性的重要性：补充材料中的一个关键负对照实验表明，如果人为地截断GPCR的一个重要胞内环（ICL3），MDPath分析出的路径就会变得模糊不清，甚至出现矛盾的信号（如在激活模拟中出现失活路径）。这证明了使用完整的、高质量的蛋白质模型进行MD模拟是获得可靠变构路径的前提。变构调节剂的影响：补充材料（图S7）还探究了钠离子和胆固醇等变构调节剂对通讯路径的影响。结果显示，这些调节剂的加入虽然会改变某些路径的权重（如增强了钠离子结合位点周围的信号），但核心的通讯通路模式保持不变，显示了变构网络的稳健性。在激酶靶点中的应用：图4：(A) ABL激酶与波舒替尼（紫色路径）和阿西米尼（蓝色路径）结合的完整视图。(B) DFG基序被变构路径稳定在DFG-out构象。(C) 远端T212残基作为正构路径的终点。图中紫色路径起始于正构ATP结合口袋（波舒替尼结合位点），蓝色路径起始于变构肉豆蔻酰口袋（阿西米尼结合位点）。两条路径都汇聚到自抑制性SH3结构域，但通过不同的机制。子图(B)显示蓝色变构路径如何稳定DFG基序（橙色）保持DFG-out构象，为阿西米尼的变构抑制机制提供分子基础。子图(C)展示远端T212残基（橙色）作为紫色正构路径的终点，解释了该位点突变如何影响ATP结合口袋抑制剂的活性。为了证明方法的普适性，作者将其应用于著名的ABL1激酶。MDPath成功识别出由正构抑制剂（波舒替尼）和变构抑制剂（阿西米尼）引发的两条截然不同的路径，并首次从动力学网络角度揭示了阿西米尼的变构抑制机制。 Q&A Q1: 这个工具对于药物研发的实际价值体现在哪里？ A1: MDPath的价值主要体现在以下几个方面：理解药物作用机制：通过可视化不同药物（如激动剂vs拮抗剂）引发的特异性通讯路径，可以深入理解其产生不同药理效应的分子基础。指导理性药物设计：识别出的路径上的关键“枢纽”残基，可以作为新的药物设计靶点，或者用于指导对现有分子的结构优化。解释耐药性突变：MDPath可以找到连接药物结合位点与远处突变位点的变构路径，从而解释为什么一个远端的突变会影响药物的疗效。发现新的变构口袋：通过分析整个蛋白的通讯网络，有可能识别出此前未被发现的、对蛋白功能至关重要的“热点”区域，这些区域可能成为全新的变构药物靶点。 Q2: MDPath的分析依赖于MD模拟，那么模拟的时长和质量对结果有什么影响？ A2: 这是一个非常关键的实际问题。模拟的时长决定了构象采样的充分性。本文使用了200 ns的模拟，这对于捕捉局部、快速的二面角运动是足够的，可以很好地分析处于一个稳定状态的通讯网络。但如果想要研究从非活性态到活性态的完整转变过程，这种慢过程就需要更长的模拟或结合增强采样方法。模拟的质量，如力场的准确性、体系构建的合理性，直接决定了轨迹的物理真实性。如果模拟本身不准确（如本文补充材料中ICL3截断的例子），那么从中分析出的任何“路径”都将是不可信的。因此，高质量、充分采样的MD模拟是MDPath分析成功的基石。 Q3: 论文中提到了对路径进行“层次聚类”，这一步的目的是什么？ A3: Dijkstra算法会找到成百上千条独立的“最优”路径。许多路径在空间上可能是高度重叠、非常相似的，它们实际上代表了同一条宏观的通讯“干道”。层次聚类的目的就是将这些相似的路径自动地分组归类。MDPath通过计算不同路径上残基原子坐标的距离来衡量路径的相似性，然后将相似的路径聚成一类。这样做的好处是，可以从纷繁复杂的数百条路径中，提炼出几条（如3-5条）最具代表性的、结构上不同的核心通讯通路（path clusters），如图3A中显示的红色和蓝色两条截然不同的路径。这极大地简化了结果的分析和可视化，让研究者能更容易地抓住主要的变构机制。关键结论与批判性总结核心结论本文成功开发并开源了一款名为MDPath的Python工具包，用于从MD模拟轨迹中系统性地识别、分析和可视化蛋白质的变构通讯路径。该方法以残基主链二面角的归一化互信息（NMI）为核心，结合图论算法，能够有效捕捉残基间的动态协同运动，并绘制出信息传递的最优路径。在GPCRs和ABL1激酶等多个重要药物靶点上的测试表明，MDPath不仅能准确识别已知的保守变构基序和激活机制，还能揭示配体特异性的信号通路。 MDPath的分析结果与实验突变数据高度吻合，能够为突变如何影响蛋白质功能提供合理的动力学机理的解释。潜在影响为药物研发领域的研究者提供了一个易于使用且功能强大的开源工具，有助于加深对药物作用机制的理解，并指导基于结构的理性药物设计。其“配体特异性”路径分析功能，为研究GPCR功能选择性、偏向性激动等前沿问题提供了新的计算视角。存在的局限性该方法目前仅考虑了主链二面角的信息，忽略了侧链运动和水分子等其他可能参与变构通讯的因素。路径识别的准确性依赖于MD模拟的充分采样。对于涉及大的构象变化的慢过程，可能需要更长的模拟或结合增强采样方法。路径的可视化和解读在一定程度上仍需要研究者的专业知识和判断。未来研究方向将侧链构象、水分子网络等更多维度的信息整合到NMI计算中，以构建更全面的通讯网络模型。将MDPath与马尔可夫状态模型（MSM）等方法结合，分析不同构象状态之间的转变路径。利用MDPath分析更大规模的MD数据库（如GPCRmd），进行高通量的变构机制探索。

Molecular Dynamics · 2025-10-08

跨越毫秒到秒级鸿沟：加权系综模拟如何捕捉”看不见”的生物动力学本文信息标题：加权系综模拟：方法、软件与应用的进展作者：Lillian T. Chong, Daniel M. Zuckerman 发表时间：2025年5月6日（ChemRxiv预印本）单位：匹兹堡大学（美国），俄勒冈健康与科学大学（美国）引用格式：Chong, L. T., & Zuckerman, D. M. (2025). WEIGHTED ENSEMBLE SIMULATION: ADVANCES IN METHODS, SOFTWARE, AND APPLICATIONS. ChemRxiv. https://doi.org/10.26434/chemrxiv-2025-jtppp 相关软件：本文主要讨论了基于 WESTPA 软件包的进展，并提及了其他实现如 wepy。摘要二十多年来，加权系综（Weighted Ensemble, WE）路径采样策略以远低于传统模拟的计算成本，实现了对罕见事件（或称跨能垒过程）路径的模拟，同时保持了严谨的动力学信息。本综述重点介绍了WE在方法和软件方面的最新进展，包括用于路径系综机理分析和高效速率估算的工具。我们展示了加权系综在一系列广泛的凝聚相过程中的成功应用，例如，微秒时间尺度的化学反应的混合量子力学/分子力学（QM/MM）模拟，以及毫秒到秒时间尺度的更慢过程的原子级模拟。这些应用涵盖了药物跨膜渗透、配体解离以及SARS-CoV-2刺突蛋白的大尺度开放等前沿领域。我们还讨论了WE策略当前面临的局限性和关键挑战，该方法尚未完全发挥其潜力。核心结论 WE是高效的罕见事件采样方法：它通过复制（分裂）和删减（合并）轨迹，能够以更低的计算成本模拟药物解离、蛋白质构象变化等低概率事件，同时严格保留动力学信息。方法学日趋成熟：近年来，WE在反应坐标优化（如机器学习辅助）、速率常数估算和不确定性量化等方面取得了显著进展，使其更加强大和可靠。软件生态系统完善：以 WESTPA 为代表的开源软件包具有高度可扩展性和互操作性，无需修改动力学引擎即可与AMBER、GROMACS、OpenMM等主流软件无缝对接，极大地促进了其应用。应用成果斐然：WE已成功应用于多个前沿领域，包括模拟秒级的SARS-CoV-2刺突蛋白开放、药物分子从深埋口袋中的解离、以及微秒级的QM/MM化学反应，揭示了实验难以企及的机理细节。背景 mindmap root((**背景与动机**)) **罕见事件采样挑战** 蛋白质折叠 **药物结合解离** 酶催化反应跨越能垒的过程 **传统MD模拟的限制** **毫秒到秒时间尺度** **计算成本高昂** 难以捕捉罕见事件 **WE策略的特色** **优胜劣汰重点培养** **动态资源分配** **无偏轨迹采样** 保留动力学信息在分子模拟的世界里，许多最引人入胜的生物学过程——如蛋白质折叠、药物分子与靶点的结合与解离、酶催化反应——都属于”罕见事件“。这意味着这些过程虽然至关重要，但在整个模拟时间尺度中，系统大部分时间都处于稳定的能量”盆地”中，而跨越能垒发生关键转变的瞬间则极为短暂和稀少。使用传统的分子动力学（MD）模拟，想要捕捉到这些事件的完整路径和动力学信息，往往需要运行长达毫秒、秒甚至更长时间的模拟，这对于目前的计算资源来说是极其昂贵甚至是不可能的。为了攻克这一难题，科学家们开发了多种增强采样和路径采样方法。其中，加权系综（Weighted Ensemble, WE）是一种尤为强大且独特的路径采样策略。与那些通过修改能量势面来加速转变的方法不同，WE的核心思想是”优胜劣汰，重点培养“。它并行地运行大量短时间的、完全标准的MD轨迹，并为每条轨迹分配一个”权重”。在固定的时间间隔后，它会评估所有轨迹的位置，智能地”克隆”那些正在向我们感兴趣的罕见区域探索的轨迹（分裂），并”删减”那些在已充分采样的区域中冗余的轨迹（合并）。通过这种方式，WE将计算资源动态地重新分配到那些”有前途”的路径上，极大地提高了采样到罕见事件的效率，同时由于每条轨迹本身是无偏的，整个过程保留了严谨的动力学信息，可以直接用来计算反应速率常数等关键物理量。经过二十多年的发展，WE方法本身、支持它的软件以及其应用范围都取得了长足的进步。关键科学问题作为一篇综述，本文旨在系统性地回答以下问题，为相关领域的研究者提供一份全面的指南和前沿展望： WE方法的核心原理是什么？它与其他路径采样方法相比有何独特的优势和固有的局限性？近年来WE方法学本身有哪些关键突破？研究者们是如何解决诸如如何定义”进展”、如何更准确地计算速率、以及如何评估结果不确定性等核心挑战的？支持WE模拟的软件生态系统发展如何？以WESTPA为代表的软件包在可扩展性、易用性和与其他主流模拟软件的兼容性方面取得了哪些进展？ WE在解决实际科学问题上取得了哪些里程碑式的应用成果？它如何帮助我们理解从病毒入侵到药物设计等一系列复杂生物过程的动力学机理？ WE方法的未来在哪里？它仍然面临哪些挑战，以及未来的发展方向将如何进一步拓展其应用边界？研究内容核心理论：加权系综（WE）模拟的”道”与”术” mindmap root((**WE核心原理**)) **基本算法** 初始化反应坐标定义箱子bins划分 **权重归一化** **动力学演化** 并行短时MD 无偏轨迹生成 **重采样操作** **分裂Splitting** 探索稀有区域克隆轨迹 **合并Merging** 删减冗余保持权重和为1 迭代循环 **动力学计算** **源-汇边界** **非平衡稳态NESS** **速率常数** **显著特点** 互操作性强算法灵活轨迹无偏连续统计严格精确 **固有局限** 物理时间尺度限制轨迹相关性问题方差挑战 WE方法的核心思想在于通过操控一个带有权重的轨迹系综，在不偏离真实动力学的前提下，高效地对罕见事件进行采样。基本算法流程：分裂与合并的艺术 WE方法通过并行管理一组带有权重的轨迹来高效探索构象空间。整个过程可以被看作一个不断迭代的循环。图1：加权系综策略示意图该图展示了一个基础的WE实现，其中构象空间被划分为固定的”箱子（bins）”，每个箱子的目标轨迹数为3。子轨迹会均分其父轨迹的权重，确保每一轮迭代中总权重为1。初始化与空间划分：首先，需要定义一个或多个”反应坐标（Progress Coordinates）“，它们是能够描述系统从初始态向目标态演化进程的变量。基于这些坐标，整个构象空间被划分成一系列离散的”箱子（bins）“。然后，从一个或多个初始构象开始，启动若干条轨迹，并为它们分配初始权重。所有轨迹的权重总和必须恒为1，即： [\sum_{i} w_i(t) = 1] 动力学演化（Evolve）：在一个迭代步中，所有轨迹都独立、无偏地进行一小段固定时间（$\tau$）的MD模拟。这个步骤是完全并行的，因此WE具有极好的可扩展性。重采样（Resampling）：这是WE的灵魂所在。在 $\tau$ 时间后，暂停所有轨迹，并根据它们所处的”箱子”进行分裂（Splitting）和合并（Merging）操作：分裂（Splitting）：当一条轨迹进入了一个很少被访问或完全空的箱子时，表明它正在探索新的、重要的区域。此时，系统会将其”克隆”成两条或多条（例如2条）子轨迹。这些子轨迹完全继承父轨迹的坐标和速度，并均分其权重（例如，权重为 $w_p$ 的父轨迹分裂成两条权重各为 $w_p/2$ 的子轨迹）。这相当于将计算资源动态地聚焦到有前途的探索路径上。合并（Merging）：当一个箱子里的轨迹数量超过了预设的目标值时，说明该区域已被过度采样，存在冗余计算。此时，系统会从中选择轨迹进行合并。例如，从箱子中随机选取两条轨迹 $i$ 和 $j$，它们的权重分别为 $w_i$ 和 $w_j$。系统会根据权重以概率 $p_i = w_i / (w_i + w_j)$ 保留轨迹 $i$，或以概率 $p_j = w_j / (w_i + w_j)$ 保留轨迹 $j$。幸存的轨迹将获得两者合并后的总权重 $w_{\text{new}} = w_i + w_j$，而被淘汰的轨迹则终止。这相当于剪除冗余的计算分支，节约资源。迭代：完成重采样后，所有”幸存”和”新生”的轨迹进入下一轮的”演化-重采样”循环，周而复始，直到达到预定的总模拟时间或目标事件被充分采样。 graph TD subgraph "方向：从左到右" direction LR A("1.初始化 一组带权重的轨迹") --> B["2.动力学演化 所有轨迹独立运行一小段时间 τ"]; B --> C{"3.重采样 (根据轨迹位置)"}; C -- "进入稀有区域" --> D["分裂 (复制有前途的轨迹)"]; C -- "进入拥挤区域" --> E["合并 (删除冗余的轨迹)"]; D --> F["进入下一轮迭代"]; E --> F; end 动力学性质计算：速率常数 WE的一个核心优势是能够直接计算动力学速率常数。这通常通过设置”源-汇（source-sink）“边界条件来实现：当一条轨迹到达我们定义的目标态（汇），它不会终止，而是被”传送”回初始态（源）并继续模拟。经过一段时间的模拟，系统会达到一个非平衡稳态（Non-Equilibrium Steady State, NESS），此时单位时间内从源到达汇的概率通量（Flux）将趋于一个稳定值，这个值就是我们要求的速率常数 $k_{AB}$。 [k_{AB} = \text{Flux}(A \rightarrow B \text{NESS})] 公式的通俗解释这个公式是WE计算速率的核心。 $k_{AB}$：是从状态A到状态B的速率常数，单位是时间的倒数（如 $\mathrm{s}^{-1}$）。 $\text{Flux}(A \rightarrow B)$：指的是单位时间内，从初始态A区域”流向”目标态B区域的总概率。在WE中，这就是所有首次到达目标态B的轨迹的权重之和除以时间间隔 $\tau$。 NESS：表示这个计算必须在系统达到非平衡稳态后进行。如图2所示，模拟刚开始时，通量会逐渐增加（瞬态），只有当进入和离开各个区域的概率流达到一种动态平衡时，测得的通量才是稳定且准确的。图2：从WE模拟流入目标态的通量估计速率常数模拟开始后，流入目标态的通量会经历一个瞬态增长期，最终达到一个平台期，即非平衡稳态，此时的通量值即为速率常数 $k$。 WE方法的显著特点与固有局限优点互操作性强：WE算法只要求能启停轨迹，因此无需修改任何MD引擎的底层代码，可以与AMBER、GROMACS、OpenMM等任何模拟软件无缝协作。这种设计使得研究者可以继续使用最熟悉、最适合其体系的MD引擎，而不必为了使用WE而去学习一个全新的模拟软件。算法灵活：WE的分箱策略、资源分配等都可以在模拟过程中动态调整，甚至可以完全抛弃”箱子”概念，而是基于轨迹间的相似度进行重采样（如REVO方案）。这种高度的灵活性使得WE能够适应各种复杂的生物分子体系。轨迹无偏且连续：WE不施加任何偏置力，每条轨迹片段都是真实的动力学路径，最终可以拼接成完整的、可用于各种机理分析的连续轨迹。这种无偏性是WE与其他增强采样方法的根本区别。统计上严格精确：理论上，WE的系综平均结果与大量传统MD模拟的结果是完全一致的。这种统计上的严谨性使得WE计算得到的速率常数等动力学可观测量具有理论上的精确性。高效并行性：WE具有极好的可扩展性，能够在数千个CPU/GPU核心上高效并行，其优秀的任务管理器和通信层设计使其能够驾驭超算级别的计算资源。计算成本显著降低：相比传统MD模拟，WE能够以远低于传统模拟的计算成本实现对罕见事件（或称跨能垒过程）路径的模拟，同时保持严谨的动力学信息。局限性 2.3 WE的内在局限性 WE方法的主要局限性源于系统固有的物理时间尺度，因此这是任何模拟真实连续轨迹系综的方法都面临的挑战。具体而言，任何感兴趣的转变过程都可以用平均过渡路径时间（average transition path time） $\langle t_{\text{TP}} \rangle$ 来表征。因此，包含 $n \gg 1$ 条轨迹的系综所需的总时间为 $n \cdot \langle t_{\text{TP}} \rangle$，这代表了在能够完全独立生成正确分布的过渡轨迹（这实际上是不可能的）的理想情况下的绝对最小计算成本。实际上，还存在一个额外的低效因子 $m > 1$（很可能 $m \gg 1$），它代表了生成独立轨迹的开销成本。因此，系综的总成本为 $m \cdot n \cdot \langle t_{\text{TP}} \rangle$，这还没有考虑 $t_{\text{TP}}$ 值可能存在的非高斯大幅度涨落。即使对于 $\langle t_{\text{TP}} \rangle \sim 10 \mathrm{ns}$ 的转变过程，在WE或其他生成连续路径系综的方法中也可能需要数微秒的轨迹数据。对 $\langle t_{\text{TP}} \rangle$ 的估计各不相同：小蛋白折叠（微秒到毫秒时间尺度）约为1-100 ns，扩散控制的蛋白-蛋白结合（微秒时间尺度）约为5 ns，蛋白-配体解离（秒时间尺度）约为100 ns。为什么高度相关轨迹会导致WE估计的可观测量（如速率常数）在不同运行之间存在高方差？统计独立性缺失：在WE中，分裂操作产生的子轨迹共享相同的历史，导致它们高度相关。这些相关轨迹不提供独立的统计信息，相当于减少了有效样本量。当多个相关轨迹贡献到同一统计量时，它们不能像独立轨迹那样有效降低方差，导致估计的不确定性增加。路径空间采样不均衡：相关轨迹倾向于探索相似的路径空间区域，使得某些重要但罕见的路径可能被低估，而常见路径则被过度采样。这种采样不均衡性会导致不同WE运行之间对同一物理量的估计出现较大波动。权重分布偏差：由于合并操作基于权重进行随机选择，高度相关的轨迹可能导致权重分布出现偏差。这种权重偏差会进一步放大估计量的方差，尤其是在长时间模拟中。收敛速度降低：相关轨迹减慢统计收敛速度，因为系统需要更长时间探索不同的路径空间。在有限的计算资源下，这可能导致不同运行之间结果差异显著。低效因子 $m$ 正是反映了这种基于相关性的低效率。在WE中，这种相关性源自基本的分裂操作。同一父轨迹的子轨迹在分裂点之前共享相同的历史，使它们高度相关。高度相关轨迹的实际后果是WE估计的可观测量（如速率常数）在不同运行之间可能存在高方差（图2）。这种基于相关性的方差和低效率可以在一定程度上得到改善，下文将详细讨论。我们还注意到，相关性使得不确定性量化更具挑战性，这也将在下文讨论。总体而言，虽然WE是一种强大而严格的方法，但并不保证在每个系统上都能很好地工作。例如，高电荷配体从蛋白受体解离是一个特别具有挑战性的压力测试；相比之下，更容易处理的应用涉及不带电配体的解离（见第5.2节）。基于系统物理性质的固有成本是显著的，这不仅对WE如此，对任何提供真实过渡路径系综的方法都是如此，即使是粗粒化模型也是如此。基于相关性的低效率也是路径采样方法的固有特性。 WE方法学的最新进展图3：WE方法学中的挑战与解决方案（a）WE模拟面临着寻找好的反应坐标、速率估计方差大和不确定性量化等挑战。（b）针对这些挑战，研究者开发了机器学习、方差最小化分箱和贝叶斯分析等解决方案。（c）这些优化方法通常需要初步模拟数据，通过分析或机器学习来指导后续的优化模拟或直接计算可观测量。近年来，研究者们从多个角度对WE方法进行了优化，主要分为两大类：优化模拟过程：反应坐标与分箱策略：这是WE实践中最关键的一环。除了依赖化学直觉，多种自动化策略被开发出来。例如，最小自适应分箱（MAB）方案能自动识别路径上的瓶颈区域并增加采样；REVO 方案则完全抛弃箱子，基于轨迹相似度进行重采样。机器学习也被用于从业已产生的高维轨迹数据中自动学习出最优的低维反应坐标，例如使用卷积变分自编码器来压缩轨迹信息。更有甚者，可以直接以最小化速率常数估计的方差为目标来优化分箱策略。优化数据分析：速率常数估算：为了解决模拟时间不足以达到稳态的问题，研究者开发了历史增强马尔可夫状态模型（haMSM），它可以从非稳态的瞬态数据中外推出稳态的速率常数。机理量化：如何从大量的路径中提取并量化”反应机理”是一个开放性问题。目前已有如LPATH等工具被开发用于对路径进行聚类和分析，以识别不同的反应通道。不确定性量化（UQ）：由于轨迹相关性，简单的统计方法不适用。目前的主流做法是进行多次独立的WE模拟，然后分析多次模拟结果之间的差异，有时还会借助贝叶斯分析来处理方差较大的情况。 WE软件的进展：以WESTPA为例 mindmap root((**WESTPA软件生态**)) **可扩展性** 数千CPU和GPU核心超算级别支持优秀任务管理器通信层设计 **互操作性** **与引擎解耦** 命令行调用 **支持主流软件** AMBER GROMACS OpenMM 无需代码修改 **数据管理** **WESTPA2.0改进** **HDF5格式** 高效存储便利重启分析数据共享优化 **未来发展** Dask任务分发减少延迟容错能力云计算支持 WESTPA（The Weighted Ensemble Simulation Toolkit with Parallelization and Analysis）是目前最活跃、功能最强大的开源WE软件包之一。高度可扩展：WESTPA能够在数千个CPU/GPU核心上高效并行，其优秀的任务管理器和通信层设计使其能够驾驭超算级别的计算资源。强大的互操作性：WESTPA设计上与动力学引擎解耦，可以像”指挥官”一样通过命令行调用任何模拟软件（如AMBER、GROMACS、OpenMM）或分析工具（如MDAnalysis、MDTraj），无需任何代码修改。数据管理优化：最新的WESTPA 2.0版本改进了数据存储框架，使用高效的HDF5格式来管理数千万个轨迹文件，极大地便利了模拟重启、数据共享和后分析。未来发展：未来的WESTPA将集成更先进的任务分发框架（如Dask），以减少延迟、增强容错能力，并更好地支持云计算平台。 WE应用的亮点成果 mindmap root((**WE应用领域**)) **病毒学** **SARSCOV2刺突蛋白** 秒级时间尺度百万原子体系 **戈登贝尔奖** 聚糖门控机制实验验证 **药物发现** **配体解离** 秒级过程不带电配体 **HIF2α靶点** 两条解离路径 **隐蔽口袋探索** 不可成药靶点药物设计新路线 **跨膜渗透** 虚拟生物利用度 **ADMETOX评估** POPC脂双层与实验一致机理洞察 **化学反应** **QMMM模拟** 微秒级反应 click反应颠覆扩散控制假设限速步骤分析 **蛋白质相互作用** 结合动力学 kon速率计算关键残基识别蛋白质折叠图4：近期WE在微秒至秒时间尺度上的应用（a）微秒级：化学反应的QM/MM模拟。（b）毫秒级：药物分子的跨膜渗透。（c）秒级：配体从深埋的受体口袋中解离。（d）秒级（百万原子体系）：SARS-CoV-2刺突蛋白的开放过程。病毒学：SARS-CoV-2刺突蛋白开放迄今为止最雄心勃勃的WE应用是对包含近百万个原子的SARS-CoV-2刺突蛋白（S蛋白）开放过程的模拟，这是一个秒级时间尺度的事件。这项工作荣获了2020年戈登·贝尔COVID-19研究特别奖。模拟不仅捕捉到了S蛋白从”关闭”到”开放”状态的完整路径，还揭示了一个前所未知的机理：位于N288位点的一个聚糖扮演了”构象门”的角色，控制着蛋白的开放。这一发现随后得到了实验的验证，包括生物层干涉测量实验和冷冻电镜（采用ManifoldEM方法生成S蛋白的大尺度运动，发现与模拟一致）。药物发现：配体解离与”隐蔽口袋”探索药物的疗效与其在靶点上的停留时间（与解离速率成反比）密切相关。WE已被成功用于模拟药物分子从靶点蛋白解离的秒级过程，迄今限于不带电配体。在一项针对癌症靶点HIF-2α PAS-B结构域的研究中，WE模拟在不知道任何先验信息的情况下，成功捕捉到了一个药物样不带电配体从其深埋的内部口袋中逃逸的两条不同路径。这些路径是以盲目方式生成的，无需任何关于解离过程的先验知识。模拟发现的构象门控残基也得到了NMR动力学实验的证实。此外，WE还能采样到在实验结构中不可见的”隐蔽口袋”，为”不可成药”靶点提供了潜在的可行药物设计路线。药物跨膜渗透：虚拟生物利用度分析 WE被用于开发预测药物被动跨膜渗透性的”虚拟实验”，这是评估药物吸收、分布、代谢、排泄和毒性（ADME/Tox）的关键性质。作为概念验证，WE模拟评估了一系列不同大小、形状和柔性的药物样胺类化合物通过模型POPC脂双层的渗透性。结果产生的渗透系数与MDCK-LE细胞系和平行人工膜渗透实验（PAMPA）的实验值一致，同时提供了转运过程的机理洞察。值得注意的是，尽管使用了被其他方法认为次优的反应坐标（膜中的z位置），WE仍成功生成了路径和速率估计，计算成本比传统MD低几个数量级。因此，WE策略对反应坐标选择的敏感性远低于基于自由能的方法。化学反应：QM/MM模拟揭示反应机理通过与混合量子力学/分子力学（QM/MM）方法结合，WE首次被用于模拟溶液中的化学反应并计算速率。在一项对叠氮化物”click反应”的研究中（叠氮阴离子与三苯甲基阳离子在乙腈-水溶液中反应），WE-QM/MM模拟不仅重现了实验速率，还颠覆了之前的”扩散控制”假设，指出反应的限速步骤是离子对中间体重排为产物的活化过程。研究还揭示了叠氮离子在阳离子苯环间增加的”爬行”与更慢的反应速率相关，这项工作突显了WE在使用混合QM/MM模型进行路径采样和动力学分析以获得更深入机理洞察方面的威力。蛋白质-蛋白质相互作用：结合动力学研究 WE已被用于研究蛋白质-蛋白质结合路径和速率常数计算，采用完全连续的显式溶剂模拟。通过模拟结合路径和解离过程，WE能够揭示相互作用界面的关键残基和构象变化。例如，WE已被用于计算基础 $k_{\text{on}}$（直接模拟柔性分子模型的蛋白-蛋白结合），以及比较无序肽及其精确预组织类似物的结合动力学。蛋白质折叠：超快折叠蛋白研究 WE已被成功应用于研究蛋白质折叠动力学和机制。例如，在对超快折叠蛋白NTL9的研究中，WE模拟揭示了改变骨架组成对折叠动力学和机制的影响。这些应用展示了WE在解决从微秒到秒原子级折叠时间的计算估计方面的能力。多尺度过程与未来展望 WE方法的应用范围正在不断扩展。除了上述应用，WE还被应用于肽跨膜渗透、脂质相分离热力学、以及大规模生物分子复合物的动力学研究。随着计算能力的提升和方法的持续改进，WE有望在更复杂的细胞环境（如呼吸道气溶胶、细菌或人类细胞质）中模拟生物分子的行为。 Q&A Q1：加权系综（WE）和其他增强采样方法（如元动力学、伞形采样）的根本区别是什么？ A1：根本区别在于是否改变系统的哈密顿量（即能量势面）。元动力学、伞形采样等方法属于偏置势（Biasing Potential）方法。它们通过在构象空间中添加一个外部的、人为的偏置势能来”填平”能量势垒，从而迫使系统更快地在不同状态间转换。这些方法能高效地计算自由能曲线，但其产生的轨迹不是真实的动力学路径，因此不能直接用来计算速率常数或分析动力学机理。加权系综（WE）则是一种路径采样（Path Sampling）方法。它不施加任何偏置力，系统在每一步都遵循自然的动力学演化。它的加速效果来自于在路径空间中对轨迹进行智能的复制和删减，即把计算资源集中到更有可能发生转变的路径上。因此，WE产生的轨迹是物理上真实的、无偏的连续路径，既可以用来计算自由能，也可以直接用来分析动力学机理和计算速率常数。 Q2：什么是好的”反应坐标（progress coordinate）”，为什么它对WE模拟如此重要？ A2：一个好的”反应坐标”是一个或一组能够有效区分反应物、产物以及过渡态的低维变量。它应该能够捕捉到系统从初始态向目标态的”进展程度”。在WE模拟中，反应坐标直接决定了”箱子（bins）”的划分，从而控制着轨迹的分裂与合并策略。一个好的反应坐标能让WE算法准确地识别出哪些轨迹正在接近反应的”瓶颈”区域（即能垒顶部），并及时在这些关键区域增加采样（分裂轨迹），从而大大提高模拟效率。相反，如果选择了一个与反应真实路径无关的坐标，WE可能会在不相关的区域浪费大量计算资源，导致收敛缓慢甚至失败。 Q3：WESTPA软件的一大亮点是”互操作性（interoperability）”，这具体指什么，为什么它很重要？ A3：互操作性指的是WESTPA能够与几乎任何现有的动力学模拟软件（如AMBER、GROMACS、OpenMM）或分析工具无缝协作，而无需对这些软件进行任何代码修改。WESTPA就像一个”总指挥”，它通过标准的命令行接口来启动、监控和停止由其他软件执行的短时间模拟任务，然后在每个迭代周期结束后收集结果并执行重采样。这一点至关重要，因为它极大地降低了使用WE方法的门槛。研究者可以继续使用他们最熟悉、最适合其体系的MD引擎，而不必为了使用WE而去学习一个全新的、功能可能不全的模拟软件。这种模块化的设计也使得更换动力学引擎或升级版本变得非常简单。关键结论与批判性总结潜在影响解锁长时程动力学：WE及其相关软件的发展，使得在原子级别上直接模拟并分析毫秒至秒级甚至更长时间尺度的生物过程成为可能，为理解药物停留时间、病毒入侵机理等关键问题提供了前所未有的工具。连接理论与实验：WE能够直接计算速率常数等动力学可观测量，这为力场的动力学性质验证提供了黄金标准，有助于推动下一代更精确的分子力场的开发。推动多尺度模拟：WE的灵活性使其不仅限于分子模拟，还可以应用于系统生物学、天气预报等更宏观的尺度，展现了其作为一种通用罕见事件采样方法的巨大潜力。研究局限性方法仍在发展中：尽管取得了巨大成功，但WE方法仍处于活跃的发展阶段。如何系统性地选择最优反应坐标、如何更精确地进行不确定性量化等问题仍是当前研究的热点和挑战。对特定体系的挑战：对于某些体系，如高电荷配体的解离，WE模拟仍然面临巨大挑战，结果的方差可能非常大，难以收敛。资源需求依然可观：虽然WE相比传统MD效率极高，但模拟秒级过程仍然需要巨大的计算资源（如SARS-CoV-2的研究），这限制了其在普通实验室的广泛应用。未来方向 QM/MM与WE的深度融合：进一步推动WE在QM/MM模拟中的应用，有望在更长的时间尺度（多微秒级）上研究酶催化和溶液中的化学反应。超长时程模拟：随着计算能力的提升和算法的持续优化，WE有望挑战秒级以上的生物过程，为研究治疗性相关的动力学事件提供更精确的速率估计。与实验数据的整合：将WE产生的路径系综与单分子实验（如FRET）或时间分辨结构生物学数据相结合，以更全面的视角揭示生物大分子的功能机理。向更复杂环境迈进：随着细胞环境的结构数据日益丰富，未来的WE模拟将不再局限于孤立的生物分子，而是能够模拟其在呼吸道气溶胶、细菌乃至人类细胞质等更真实、更拥挤环境中的行为。

Molecular Dynamics

Mendelevium

Contact

Sampling & Analysis